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INTRODUCTION 


The  Navy  operates  and  maintains  many  shore-based  communication 
facilities  that  employ  large  insulators  in  their  antenna  arrays.  Due  to 
various  problems  with  these  insulators,  an  investigation  of  insulator 
failures,  alternative  insulator  design,  and  new  concept  development  was 
begun. 

Electrical  difficulties,  such  as  arcing,  corona,  and  heating  or 
burning  of  insulators,  are  associated  with  the  presence  of  an  extremely 
high  electric  field.  Therefore,  when  evaluating  old  and  new  insulator 
designs,  it  is  helpful  to  know  the  electric  field  distribution  asso- 
ciated with  each  insulator.  The  geometric  complexity  of  insulators  and 
their  associated  electric  fields  precludes  the  use  of  classical  analytical 
methods  to  determine  these  fields.  This  report  describes  a computer 
program  developed  to  calculate  and  plot  the  electric  potential  distribu- 
tion of  insulator  configurations. 

Many  antenna  insulators  have  geometries  that  are  axially  symmetric; 
that  is,  the  insulator  appears  the  same  regardless  of  how  it  is  turned 
about  its  axis  of  symmetry.  The  same  is  true  of  the  electric  fields 
associated  with  such  an  insulator.  Thus,  the  problem  can  be  reduced 
from  one  involving  three  dimensions  to  one  of  two  dimensions. 

The  program  calculates  the  potential  distribution  for  multidielec- 
tric media  by  numerically  solving  Laplace's  equation.  The  potential 
distribution  is  found  in  a direct  manner  by  solution  of  simultaneous, 
finite-difference  equations.  Computer  drawings  of  equipotential  lines 
can  then  be  made. 

The  program  is  capable  of  solving  many  other  types  of  problems, 
including  electric  flux  lines,  mechanical  stress,  and  temperature  dis- 
tribution. Geometries  that  are  strictly  two  dimensional  can  be  treated 
by  making  a minor  modification  to  the  program. 

This  report  describes  the  mathematics  behind  and  organization  of 
the  computer  program.  A detailed  users  guide  is  not  given  here,  although 
one  is  being  prepared.  A listing  of  the  computer  statements  is  given  in 
the  Appendix. 


BACKGROUND 

The  electric  fields  of  high  voltage  insulators  satisfy  Maxwell's 
equations.  But  all  of  Maxwell's  equations  do  not  need  to  be  solved  in 
order  to  obtain  the  desired  information  about  the  insulator  fields.  In 
studying  the  electrical  breakdown  of  Insulator  configurations,  it  is  the 
electric  field  intensity  which  is  of  Interest.  The  electric  field 
intensity,  E,  consists  of  two  components  as  shown  in  Equation  1. 
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The  scalar  potential  function,  V,  is  related  to  the  voltage  applied 
across  the  insulator,  and  the  voltage  gradient,  W,  at  a given  point  is 
related  to  the  dimensions  of  the  insulator  and  the  distance  from  the 
electrodes  of  the  insulator.  The  other  component  is  the  rate  of  change 
of  the  magnetic  vector  potential.  A,  which  is  a function  of  the  current 
distribution  on  the  antenna,  the  geometry  of  the  insulator  electrodes, 
and  the  distance  from  the  insulator  electrodes. 

Since  insulators  consist  of  one  or  more  dielectric  materials  that 
are  different  from  air,  it  is  necessary  to  apply  the  constituitlve 
relation  to  Equation  1 to  get  the  displacement  flux  density,  D. 


where  e is  the  dielectric  constant  of  the  material  in  which  D is  evalu- 
ated. 

In  a region  that  contains  no  free  electric  charge,  such  as  the 
region  around  the  insulator  electrodes.  Maxwell's  divergence  equation 
for  the  electric  field  is  written  as 


V*D  = 0 


(3) 


Applying  this  condition  to  Equation  2 yields 


3t; 


or 


V»D  = - e V#W  + V*  = 0 


V,D  = - e ^ V^V  + 1^  (V»A)j  = 


V»A  is  defined  by  the  Lorentz  condition.  Equation  6, 


(4) 


(5) 


VA  + + UOV  = 0 


(6) 


where  p and  a are  the  permeability  and  conductivity  of  the  medium, 
respectively.  Using  Equation  6 to  substitute  for  7»5,  Equation  5 
becomes 
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V*D  = 
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Since  V,  in  the  steady-state  operating  condition,  is  of  the  form 


V = V (r,  z)  e 
o 


jwt 


(9) 


at  any  given  point  in  space  for  some  V^(r,  z),  the  time  derivatives  of  V 


are  given  by 
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where  w is  the  angular  frequency  of  the  system.  Table  1 shows  approxi- 
mate values  for  the  constants  in  Equation  8 for  the  space  between  the 
electrodes  of  most  types  of  VLF  antenna  Insulators. 


Table  1.  Constants  for  Regions  Near  Electrodes  of 
Most  VLF  Antenna  Insulators 


Constant 


Symbol 


Approximate  Value 


Angular  frequency 

0) 

2 

X 

10^  sec  ^ 

Dielectric  constant 

e 

1 

X 

10“^^  F/m 

Permeability 

u 

1 

X 

10~^  H/m 

Conductivity 

a 

1 

X 

10  (ohm-m)  ^ 

Rewriting  Equation  8 using  Equations  10  and  11 
V»D  = - e 


2 2 
VV  + (yew  - jya  (d)V 


= 0 


(12) 
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But 


M e J - jpao)  = 4x  lO"^  - j 2 x 10“^^  = 4 x lO"^  (13) 

This  implies  that  the  jyoaiV  term  may  be  ignored  when  considering  the 
steady-state  voltage  of  most  insulators,  since  it  is  much  smaller  than 
the  term.  Equation  12  then  becomes 

7«D  = - g(V^V  + V)  = 0 (14) 


where 


2 / 10-7 

p e 0)  = 4 X 10 


For  the  purpose  of  evaluating  Equation  14  it  is  sufficient  to  approx- 
imate the  VLF  antenna  insulator  as  a parallel-plate  capacitor  shown  in 
Figure  1.  In  the  cylindrical  coordinate  system  Equation  14  is  written 
as 

Since  the  problem  has  cylindrical  symmetry,  V is  independent  of  <f>  and 
87v/3(()^  = 0.  To  facilitate  easy  analysis,  the  problem  can  be  simplified 
by  considering  the  region  near  the  middle  of  the  capacitor,  far  from  the 
edges.  In  this  region,  V is  independent  of  r,  and  9V/9r  = 0.  Equa- 
tion 14  then  becomes: 


V = 0 

9z'^ 


This  equation  has  a simple  solution  which  can  be  easily  analyzed. 


V = — : — : — V sin  k z 
sin  k d 


For  a VLF  insulator 


-4  -1 

k = 6 X 10  meter 


d = 1 meter 
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V 

o 

k d 


k 2 


(19) 


Equation  19  is  the  solution  to  Laplace's  equation: 

V^V  = 0 (20) 

2 

It  can  be  shown  that  the  k V term  in  Equation  14  and  16  has  a minimal 
effect  on  the  solution. 

Thus,  for  the  size  and  frequency  range  of  VLF  antennas  the  potential 
distribution,  V,  can  be  found  to  sufficient  accuracy  by  solving  Laplace's 
equation. 

The  electric  field  intensity  consists  of  another  component  due  to 
the  magnetic  vector  potential  as  shown  in  Equation  1.  In  analyzing  the 
effect  of  this  component,  the  region  at  the  surface  of  an  electrode, 
as  in  Figure  2,  will  be  considered.  This  is  appropriate  since  the 
electric  field  intensity  will  be  maximum  near  the  electrode,  and  it  is 
this  region  which  will  be  most  vulnerable  to  electrical  breakdown. 
Equation  14  is  integrated  over  a small  volume  element  through  which  the 
electrode  surface  passes. 


J ,.D 

vol 


' / 

vol 


(V^V  + k^V) 


dv 


(21) 


The  V»D  and  V V terms  in  Equation  21  are  converted  to  surface  integral 
terms  by  application  of  the  divergence  theorem. 


(22) 


If  the  region  is  considered  to  be  an  infinitesimal  volume  element  so 
that  D,  VV,  and  V are  invariant  throughout  the  region,  the  equation 
becomes 
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where  Dq,  VVq,  Vq  are  values  at  the  surface  of  the  electrode.  Completing 
the  integration  yields: 

5 •AS  = - e(VV  "AS  + V AS  M)  (24) 

o o o 

The  approximate  magnitudes  of  the  various  factors  in  the  right-+iand  side 
of  Equation  24  for  VLF  antennas  are  given  in  Table  2.  The  voltage 
gradient  value  is  an  average  value  and  is  conservative  when  compared  to 
the  peak  gradients,  which  are  of  greatest  interest.  Therefore, 

jvV^.ASl  = |as1  (lx  10^  V/m) 

|k^  AS  A«,l  = 1 AS  Ail  I (2  X lO"^  V/m^) 

Keeping  in  mind  that  Ail  is  an  infinitesimal  length,  it  can  be  seen  that 
the  k^  Vq  term  is  insignificant  when  compared  to  the  voltage  gradient 
term.  Equation  24  can  be  approximated  quite  accurately  by 

5 ‘AS  = - e VV»  AS  (25) 

o o 

or,  simplifying 

D 

E = -2  - vv  (26) 

o e o 


Table  2.  Factors  in  an  Infinitesimal  Volume  Element  at  an 
Electrode  of  a VLF  Antenna  Insulator 


Factor 

Symbol 

Approximate  Value 

Voltage  gradient 

VV 

o 

1 X 10^  V/m 

Voltage 

V 

5 X 10^  V 

o 

-4  -1 

Wave  number 

k 

6 X 10  m 

In  summarizing,  the  electric  field  intensity  near  the  electrodes  of 
an  insulator  is  very  nearly  equal  to  the  voltage  gradient  at  that  point; 
the  component  contributed  by  the  magnetic  vector  potential  in  Equation  1 
can  be  neglected.  The  electric  field  intensity  can  be  interpreted  from 
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the  potential  distribution  by  Equation  26.  For  example,  if  equipotential 
lines  are  plotted,  the  regions  where  the  lines  are  most  dense  are  the 
regions  of  highest  potential  gradient  and  highest  electric  field  inten- 
sity. It  is  easier  to  visualize  the  electric  fields  of  Insulators  in 
this  manner  than  by  plotting  electric  field  lines  themselves.  For  this 
reason,  the  computer  program  is  designed  to  calculate  and  plot  the 
potential  distribution. 

COMPUTER  SOLUTION  METHOD 
Potential  Distribution  Problem 

As  described  above,  the  potential  distribution  near  an  insulator 
can  be  found  accurately  by  solving  Laplace's  equation  (Equation  20). 
Laplace's  equation  can  be  rewritten  in  the  form  of  Equation  14,  ignoring 
the  k'^V  term,  since  its  effect  on  the  solution  is  negligible: 

e(V^V)  = v«D  = 0 (27) 

Coordinate  Grid  System.  Finite-difference  techniques  are  used  to 
solve  the  equation  numerically.  The  geometry  of  the  region  near  the 
antenna  insulator  is  described  according  to  a grid  system,  such  as  the 
! one  shown  in  Figure  3.  The  grid  is  drawn  in  a plane  passing  through  the 

[ axis  of  symmetry  of  the  insulator  so  that  a section  view  of  the  insulator 

; can  be  defined  in  terms  of  the  grid  points.  The  boundary  of  a material 

I (electrode  or  dielectric)  is  approximated  by  a series  of  line  segments 

f that  connect  adjacent  grid  intersection  points.  These  line  segments  can 

; be  vertical,  horizontal,  or  diagonal.  All  line  segments  must  start  and 

• end  on  grid  intersection  points,  and  diagonals  must  connect  two  points 

J that  are  formed  by  the  intersection  of  adjacent  horizontal  and  vertical 

grid  lines.  By  properly  choosing  the  spacing  of  grid  lines  and  using 
combinations  of  horizontal,  vertical,  and  diagonal  line  segments,  the 
t material  boundaries  of  an  insulator  can  be  described  quite  accurately. 

' The  ability  to  use  variable  grid  density  is  an  important  feature. 

If  the  entire  grid  region  were  composed  of  equally  spaced  grid  lines, 

I the  total  number  of  grid  lines  would  be  determined  by  the  density 

required  to  accurately  describe  the  smallest  feature  of  the  Insulator 
I configuration.  This  often  results  in  an  excessively  large  number  of 

grid  lines  in  regions  where  they  are  not  necessary.  Allowing  a variable 
grid  density  makes  it  possible  to  use  a high  grid  density  only  in  regions 
1 where  it  is  necessary  to  achieve  an  accurate  description  of  the  Insulator 

or  where  detailed  information  on  the  potential  distribution  solution  is 
required.  The  result  is  a substantial  savings  in  computer  memory  require- 
f ments  and  program  execution  time.  Also,  if  a uniform  grid  size  is  used, 

it  is  often  necessary  to  compute  two  or  more  solutions,  using  a finer 
I grid  with  each  new  solution  in  order  to  obtain  the  desired  accuracy  in  a 

I given  region  of  interest.  When  using  a variable  grid  density,  one 

j solution  is  usually  all  that  is  required. 
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Formulation  of  Numerical  Equations.  Equation  27  is  integrated  over 
a volume  element  surrounding  each  grid  intersection  as  shown  in  Figure  4. 


e V V dv 


■ /’•” 


dv  = 


A surface  integral  is  obtained  by  applying  the  divergence  theorem  to 
Equation  28. 


^ V«D  dv  = ^ D< 


ds  = ® e VV»ds  = 0 


This  is  an  integral  of  the  displacement  flux  densities  crossing  the 
surface  areas  of  the  volume  element.  Equation  29  is  applied  to  the 
surfaces  of  the  volume  element  in  Figure  4.  Due  to  cylindrical  symmetry, 
the  voltage  gradient  in  the  (|)  direction  is  zero.  Thus,  the  integrals 
over  the  vertical  faces  in  the  (p-plane  are  also  zero.  The  gradients 
across  the  remaining  four  faces  are  generally  not  zero.  These  remaining 
faces  are  chosen  such  that  they  pass  through  the  midpoints  of  the  line 
segments  joining  the  central  point,  whose  potential  is  (t)Q,  and  the 
corresponding  horizontally  or  vertically  adjacent  point,  (l>j  through  <p^. 
Each  of  the  four  faces  can  be  divided  by  a boundary  which  separates  two 
different  dielectric  materials.  There  are,  then,  a total  of  eight  possi- 
ble dielectrics  and  eight  different  homogeneous  segments  through  which 
the  flux,  qj^  if  may  pass. 

In  calculating  the  integral  in  Equation  29  numerically,  the  flux 
crossing  the  surfaces  is  approximated  by  a finite  difference  value  for 
the  voltage  gradient.  The  normal  component  of  flux  passing  through  a 
given  face  is  assumed  to  be  uniform  across  that  face.  For  example. 


Since  this  normal  component  of  the  voltage  gradient  is  perpendicular  to 
the  surface  through  which  it  passes,  the  dot  product  in  Equation  29  is 
reduced  to  simple  multiplication.  The  integral  in  Equation  29,  then,  is 
the  sum  of  eight  terms,  each  of  which  is  a product  of  a displacement 
flux  density  and  the  area  through  which  the  flux  passes.  The  form  in 
which  this  equation  is  used  in  the  program,  with  the  coefficients  of  the 
((i ' s grouped  together  and  a factor  of  it  dropped,  is 
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0 (31) 


Organization  of  Equations.  There  is  one  equation  of  the  form  of 
Equation  31  for  each  grid  intersection  point.  If  the  grid  consists  of  m 
horizontal  and  n vertical  lines,  as  in  Figure  3,  there  will  be  mn  = m x n 
equations.  These  simultaneous  equations  can  be  written  in  matrix  form 
as 


A $ = B 


(32) 


where  A is  an  mn  by  mn  constant  matrix,  $ is  a vector  of  mn  unknowns, 
and  B is  an  mn-vector  of  constants.  By  carefully  choosing  the  order  of 
the  mn  equations,  A can  be  made  a banded,  symmetric  matrix  as  shown 
schematically  in  Figure  5.  The  ordering  of  the  equations  is  indicated 
by  the  subscripts  of  the  (}>'s  in  Figure  3.  The  r*-^  equation,  whose 
coefficients  appear  in  the  r^^  row  of  the  A matrix,  involve  at  most  five 
grid  points  whose  subscripts  fall  between  r - n and  r + n,  inclusive. 
This  means  that  the  bandwidth  of  A is  never  larger  than  2n  + 1,  with  n 
elements  on  each  side  of  the  diagonal.  Furthermore,  this  arrangement 
scheme  yields  a symmetric  matrix  so  that  only  the  diagonal  elements  and 
the  elements  on  one  side  of  the  diagonal  need  to  be  stored  in  the  com- 
puter in  order  to  know  the  entire  matrix.  That  means  storage  for  only 
mn  X (n  + 1)  elements  is  required  instead  of  mn  x mn  if  all  elements 
were  needed.  This  ordered  numbering  system  results  in  a substantial 
savings  in  computer  storage.  The  bandwidth  is  narrowest  and  memory 
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requirements  s—allest  when  n is  less  than  m.  Such  is  the  case  for  most 
axially  symmetric  antenna  insulators,  which  tend  to  have  diameters  that 
are  smaller  than  their  lengths.  This  usually  results  in  the  need  for 
fewer  vertical  grid  lines  than  horizontal  grid  lines.  Consequently, 
the  computer  program  is  arranged  such  that  the  vertical  grid  lines  are 
parallel  to  the  axis  of  symmetry. 

There  are  few  nonzero  elements  in  B.  These  nonzero  elements  cor- 
respond to  equations  which  involve  grid  points  that  are  of  a known 
specified  potential.  For  example,  if  a given  (pj.  is  a point  on  an  elec- 
trode of  potential  Vq,  Equation  31  is  not  used,  and  the  equation  for  that 
point  is  simply 


(33) 


where  V becomes  an  element  of  the  B vector, 
o 

Solution  of  Equations.  Once  all  of  the  mn  equations  are  formulated, 
the  matrix  equation  (Equation  32)  is  converted  to  reduced  echelon  form 
by  the  process  of  Gauss-Jordan  elimination.  The  ((>  solutions  are  then 
calculated  by  back  substitution  into  the  reduced  set  of  equations. 

Boundary  Conditions.  There  are  several  grid  points  that  must  be 
treated  specially.  One  of  these,  in  which  a grid  intersection  is  a 
point  on  a conductor  or  electrode  of  known  potential,  has  already  been 
discussed.  Points  along  the  edge  of  the  grid  must  also  be  given  special 
consideration.  If  the  potentials  along  the  grid  boundary  are  known, 
those  points  can  be  treated  as  described  above  for  Equation  33.  If  the 
grid  boundary  potentials  are  unknown,  other  methods  must  be  used  since 
these  boundary  points  have  only  two  or  three  neighboring  points  instead 
of  four  as  in  the  general  case  of  Figure  4. 

For  a grid  boundary  that  coincides  with  the  axis  of  symmetry,  there 
are  four  grid  points  involved  as  in  Figure  6.  The  same  arguments  apply 
as  for  Figure  4,  except  that  there  are  only  four  dielectric  segments 
through  which  flux  passes.  The  displacement  flux  density  through  each 
segment  is  multiplied  by  the  surface  area  through  which  the  flux  passes. 
With  the  sum  of  these  four  terms  set  equal  to  zero,  a factor  of  v 
dropped,  and  the  coefficients  of  the  (f's  grouped  together,  an  equation 
similar  to  31  appears  that  corresponds  to  Equation  29. 
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In  many  problems  the  bottom  horizontal  grid  line  represents  a plane  of 
symmetry.  For  example,  if  the  two  ends  of  a guy  line  insulator  are 
mirror  images  of  each  other,  the  plane  which  passes  through  the  midpoint 
of  the  insulator  and  is  perpendicular  to  the  axis  of  symmetry  is  a plane 
of  symmetry.  The  method  of  Images  can  be  applied  to  simplify  the  problem 
and  obtain  the  boundary  condition.  Since  both  ends  of  the  Insulator  are 
identical,  the  solution  needs  only  to  be  calculated  from  the  plane  of 
symmetry  toward  one  end.  The  boundary  condition  at  the  plane  of  symmetry 
can  be  satisfied  by  assigning  all  points  on  the  corresponding  grid  line 
a potential  equal  to  the  value  midway  between  the  potentials  of  the  two 
insulator  ends. 

In  all  cases  where  an  explicit  boundary  potential  is  unknown, 
another  means  of  satisfying  the  boundary  condition  requirements  must  be 
found.  A simple  method  of  obtaining  the  boundary  condition  is  to  deter- 
mine an  approximation  for  a particular  electric  flux  line  that  exists  at 
a large  distance  from  the  region  of  greatest  Interest.  Figure  7 depicts 
such  a flux  line  for  a guyline  insulator  which  is  not  symmetrical  about 
its  midpoint.  The  lack  of  end-to-end  symmetry  of  the  insulator  assembly 
requires  that  three  of  the  grid  boundaries  (top,  bottom,  and  right-hand 
sides)  must  be  treated  in  this  manner.  The  location  of  the  distant  flux 
line  can  be  found  easily  by  analytical  methods,  since  most  field  problems 
degenerate  to  a very  simple  problem  when  only  the  solution  at  a distant 
point  is  required.  In  the  case  of  a guyline  insulator,  the  assembly 
appears  electrically  as  a simple  needle  gap  when  viewed  from  a large 
distance.  The  needle  gap  problem  is  one  for  which  the  fields  are  easily 
calculated.  A flux  line  is  by  definition  a line  of  direction  of  the 
electric  field.  The  flux  is  along  and  parallel  to  the  line,  and  no  flux 
crosses  the  line.  If  the  grid  intersection  labeled  (()  in  Figure  4 were 
located  on  or  inside  this  flux  line,  it  would  experience  no  displacement 
flux  density  from  outside  the  flux  line.  This  can  be  described  mathe- 
matically by  assigning  e = 0 for  the  region  outside  the  flux  line. 

There  is  no  known  material  having  a dielectric  constant  of  zero;  this  is 
merely  a mathematical  technique  that  effectively  satisfies  the  boundary 
condition.  As  a result,  the  points  beyond  the  flux  line,  in  the  region 
where  the  dielectric  constant  is  zero,  have  no  bearing  on  the  potential 
solution  in  the  interior  region,  and  it  is  not  necessary  to  know  the 
potentials  along  the  outer  grid  boundary. 

Treatment  of  Conductors.  Conductors  are  fundamentally  different 
from  dielectrics  and  must  be  given  special  treatment.  If  the  potential 
is  known,  as  with  insulator  electrodes,  the  grid  points  in  that  conductor 
are  treated  as  in  Equation  33.  But  some  insulator  assemblies  have  metal 
components  that  are  between  the  electrodes  and  are  insulated  from  any 
fixed  potential.  An  example  of  such  a "floating  potential  conductor"  is 
the  pin  that  is  used  to  connect  one  or  more  suspension  insulators  together 
to  form  a chain  of  insulators.  The  potential  of  such  a conductor  is  not 
known  and  cannot  be  specified  on  input;  it  must  be  calculated  during  the 
execution  of  the  program.  This  is  done  by  using  the  same  method  as  for 
other  points,  but  a very  high  dielectric  constant  is  used  to  describe 
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Figure  7.  Boundary  condition  at  a distant  electric  flux  line. 
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the  region  of  the  conductor.  If  the  relative  dielectric  constant  is 
nearly  infinite  when  compared  to  the  other  dielectrics  in  the  problem, 
it  is  equivalent  to  requiring  that  the  electric  field  be  zero  throughout 
that  region.  That  is,  all  points  would  be  at  the  same  potential,  as  in 
a conductor.  This  method  of  representing  floating  potential  conductors 
makes  it  possible  to  solve  the  problem  using  Laplace's  equation  instead 
of  Poisson's  equation  as  would  be  necessary  if  the  free  charge  distribu- 
tion on  the  conductors  were  considered.  The  theory  behind  the  method  is 

illustrated  in  Figure  8.  As  the  dielectric  constant  of  the  sphere  is  ■ 

increased  in  views  (a)  through  (d) , it  begins  to  appear  more  like  a 

conductor  until,  in  Figure  8d , all  points  in  the  region  are  at  the  same  ; 

potential.  The  exact  magnitude  of  the  dielectric  constant  that  is  i 

necessary  to  give  satisfactory  results  may  depend  upon  the  particular  j 

computer  being  utilized.  Using  a value  approximately  10®  times  larger  | 

than  the  other  dielectric  values  involved  in  the  problem  has  shown  to  be  j 

successful.  i 


Computer  Plotting  of  Solutions 

Calcomp  plotting  is  used  to  give  a pictorial  representation  of  the 
equipotential  lines.  From  an  equipotential  plot  it  is  easy  to  determine 
the  regions  of  highest  electric  field. 

To  plot  a given  equipotential  line,  a systematic  search  of  the  grid 
is  executed  to  determine  the  coordinates  of  points  having  the  given 
potential.  A starting  point  for  the  line  is  determined  by  successively 
searching  along  each  grid  line  until  two  points  are  found  such  that 
their  potentials  bracket  the  potential  value  to  be  plotted.  The  search 
for  this  initial  point  on  the  equipotential  line  may  include  horizontal 
or  vertical  grid  lines  or  both.  An  example  of  a vertical  grid  line  is 
shown  in  Figure  9,  where  points  A through  F represent  intersections  with 
horizontal  grid  lines.  If  a starting  point  for  the  equipotential  line 
of  value  50.0  were  desired,  it  would  be  found  in  the  interval  between 
points  C and  D.  The  location  of  the  intersection  of  the  equipotential 
line  with  this  vertical  grid  line  is  determined  approximately  by  inter- 
polation between  points  C and  D.  Several  methods  of  Interpolation  may 
be  used.  The  simplest  is  first  degree  or  linear  interpolation,  which 
assumes  a linear  variation  in  potential  between  C and  D.  A slightly 
more  sophisticated  method  involves  using  the  data  from  points  B and  E in 
addition  to  C and  D.  It  is  then  possible  to  generate  a third  degree 
polynomial  approximation  for  the  potential  distribution  between  points  B 
and  E.  This  polynomial  can  then  be  solved  to  determine  an  approximate 
vertical  coordinate  for  the  50.0  value.  The  type  of  interpolation  that 
is  employed  depends  upon  the  accuracy  required  and  the  relative  spacing 
of  the  grid  lines  with  respect  to  the  smoothness  of  the  equipotential 
curves . 


tJ-gure  9.  Section  of  vertical 
grid  line  with  potential 
solution  at  intersection 
points. 


Figure  10.  Sample  grid  rectangle  and  potential  solution  at  inter 
section  points. 


Once  an  initial  point  of  the  equlpotential  line  is  obtained,  the 
next  point  is  found  by  checking  the  remaining  sides  of  the  rectangle 
formed  by  adjacent  grid  lines  on  one  side  of  the  line  segment  CD,  as 
shown  in  Figure  10.  Since  an  equlpotential  line  must  be  continuous 
throughout  the  grid  region,  if  it  enters  a rectangle  through  one  side. 
Such  as  CD,  it  must  then  exit  through  one  of  tlie  other  three  sides. 

This  assumes  that  the  grid  network  is  fine  enough  so  that  the  equipo- 
tential  lines  cannot  enter  and  exit  on  the  same  side  of  a rectangle. 

This  requirement  is  satisfied  if  the  grid  is  designed  by  the  user  to  be 
fine  enough  in  regions  where  the  electric  field  is  sharply  curved  or 
contorted.  Subsequent  equlpotential  points  are  found  by  searching 
successive  rectangles.  For  example,  in  Figure  10  the  next  rectangle  to 
be  searched  would  be  the  rectangle  above  line  segment  DC.  If  the  desired 
potential  occurs  exactly  at  a corner  of  a rectangle,  the  potential  of 
that  grid  intersection  is  shifted  by  a finite,  but  small.  Increment, 

This  avoids  numerical  difficulties  that  would  occur  if  the  equlpotential 
line  were  to  go  through  a corner  of  a rectangle  instead  of  one  of  the 
sides.  The  incremental  shift  is  negligible  so  that  the  appearance  of 
the  plot  is  unaffected. 

Tlie  search  process  is  repeated  so  that  the  equlpotential  is  traced 
throughout  the  entire  grid  until  the  line  either  ends  at  the  grid  bound- 
ary or  circles  back  to  the  initial  point.  The  equlpotential  curves  are 
drawn  by  connecting  these  points. 

Tliere  are  several  formats  in  which  the  equlpotential  plot  can  be 
drawn.  The  equlpotential  lines  can  be  plotted  thoughout  the  entire  grid 
region,  or  small  segments  of  the  grid  can  be  enlarged  and  plotted.  The 
half-section  view  of  one  side  of  the  axis  of  symmetry  can  be  "rotated" 
about  the  axis  of  symmetry  so  that  a full  section  view  of  the  insulator 
can  be  plotted  in  one  figure.  Such  a plot,  showing  the  insulator  on 
both  the  right-  and  left-hand  sides  of  the  axis  of  symmetry  is  often 
easier  to  Interpret  than  if  only  one  side  were  plotted.  This  can  be 
taken  a step  farther  if  the  insulator  assembly  to  be  plotted  has  end-to- 
end  symmetry  so  that  the  two  ends  are  mirror  images  about  a plane  of 
symmetry.  The  view  of  one  end  may  be  "reflected"  across  the  plane  of 
symmetry  such  that  a full-length  view  of  the  insulator  is  plotted.  The 
method  used  to  interpolate  potentials  between  grid  intersections  will 
affect  the  appearance  of  the  plot  in  regions  where  grid  lines  are  sparse. 
Wliere  linear  interpolation  is  unsatisfactory,  polynomial  interpolation 
is  often  adequate  to  yield  smooth  plot  curves. 

Examples 

Fig.,i'e  11  shows  a guy  line  insulator  having  end-to-end  symmetry. 

The  potential  distribution  need  only  be  calculated  for  one-fourth  of  the 
section  view  shown  in  Figure  11.  The  solution  for  the  remaining  three 
quadrants  can  be  obtained  from  the  first  quadrant  solution  using  symmetry. 
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The  coordinate  grid  for  one  quadrant  of  the  insulator  assembly  is 
shown  in  Figure  12,  The  grid  lines  are  very  dense  near  the  origin 
(where  the  insulator  is  located),  since  this  Is  the  region  where  the 
greatest  accuracy  is  required.  The  grid  lines  are  less  dense  in  regions 
away  from  the  insulator  where  the  electric  field  is  smooth.  The  distant 
boundary  condition  is  satisfied  by  the  curved  dielectric  boundary  which 
is  an  approximation  of  an  electric  flux  line.  The  dielectric  constant 
is  unity  everywhere  inside  the  boundary  that  is  occupied  by  air.  The 
dielectric  constant  is  zero  outside  this  boundary  to  simulate  the  lack 
of  displacement  flux  density  crossing  the  flux  line.  The  boundary 
conditions  for  the  bottom  and  left-hand  sides  are  satisfied  by  the 
symmetry  of  the  configuration. 

Figure  13  is  the  equipotential  plot  for  the  quadrant  of  the  insulator 
shown  in  Figure  12.  The  voltage  between  adjacent  equipotential  lines  is 
2.5%  of  the  voltage  between  the  electrodes  of  the  insulator  in  Figure  11. 
In  general,  the  equipotential  values  that  are  to  be  plotted  can  be 
chosen  as  necessary  to  give  an  adequate  picture  of  the  potential  distri- 
bution. 

If  a more  detailed  view  of  the  region  near  the  insulator  (or  any 
other  region)  is  desired,  an  enlarged  view  can  be  obtained.  One  way  to 
do  this  is  to  plot  the  same  data  on  a larger  scale.  Figure  14  is  the 
grid  in  the  region  of  the  insulator  in  Figure  12  which  has  been  plotted 
on  a larger  scale. 

The  equipotential  plot  for  the  region  of  Figure  14  is  given  in 
Figure  15.  This  plot  was  made  using  a portion  of  the  same  data  that  was 
used  to  plot  Figure  13. 

If  the  grid  is  not  of  sufficient  density  in  the  region  of  the 
insulator,  it  may  not  be  possible  to  obtain  an  accurate  enlarged  plot 
from  that  set  of  data.  It  is  then  necessary  to  obtain  another  potential 
solution  using  a denser  grid  system  that  includes  only  the  region  to  be 
enlarged.  The  boundary  conditions  for  this  second  grid  can  be  obtained 
from  the  potential  solution  for  the  original  grid.  Not  all  of  the 
boundary  grid  points  in  the  denser  grid  will  have  existed  in  the  coarse 
grid  and  their  potentials  will  not  be  explicitly  known.  The  potentials 
of  those  boundary  grid  points  that  did  not  exist  in  the  original  grid 
can  be  found  by  interpolation  between  points  of  known  potential. 

It  is  often  desirable  to  obtain  a full  view  of  an  insulator  plot, 
rather  than  a view  of  one-half  or  one  quadrant.  The  additional  views 
can  be  plotted  by  repeating  the  plotter  pen  commands  two  or  four  times 
with  the  appropriate  inversions  to  achieve  the  full  view  of  the  desired 
region.  Figure  16  is  such  a plot  of  the  guy  line  insulator. 

The  suspension  insulator  string  in  Figure  17  is  an  example  of  an 
insulator  that  does  not  have  end-to-end  symmetry  and,  hence,  no  plane  of 
symmetry  midway  between  the  end  electrodes.  In  this  problem,  the  bound- 
ary conditions  on  three  sides  are  unknown.  By  making  an  approximation 
for  the  flux  line  shown,  these  boundary  conditions  are  satisfied. 
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Figure  13.  Equipotential  plot  for  one  quadrant  of  guy  line 
insulator. 


Figure  18  is  an  enlarged,  full-section  view,  equipotential  plot  of 
the  suspension  insulator  chain.  The  voltage  between  each  pair  of  adja- 
cent equipotential  lines  is  2%  of  the  total  voltage  across  the  end 
electrodes.  The  potentials  of  the  two  end  electrodes  may  be  arbitrarily 
specified,  but  the  potentials  of  the  interior  electrodes  are  unknown 
prior  to  execution  of  the  program.  These  conductors  must  be  treated  as 
dielectrics  having  an  infinite  dielectric  constant.  Their  potentials 
are  then  calculated  in  the  same  manner  as  for  all  other  points. 

Electric  Flux  Lines 

In  addition  to  the  equipotential  lines,  the  direction  lines  of  the 
electric  field  also  give  information  about  an  electric  field.  These 
lines  of  electric  flux,  which  are  orthogonal  to  the  equipotential  lines, 
can  also  be  found  by  solution  of  Laplace’s  Equation  35. 

V^U  = 0 (35) 


The  computer  program  will  calculate  a value  of  U for  each  grid  point  in 
a manner  similar  to  that  for  electric  potentials.  However,  unlike  the 
potential  solution,  the  numerical  values  of  U have  no  direct  signifi- 
cance. A line  of  equal  U value  is  an  electric  field  direction  line, 
but,  in  general,  there  is  no  clear  relationship  between  the  spacing  of 
these  flux  lines  and  the  values  of  U along  the  lines.  It  is  somewhat 
difficult  to  determine  the  values  of  U that  should  be  plotted  to  give  an 
accurate  picture  of  the  electric  field  intensity.  This  is  in  contrast 
to  the  plotting  of  equipotential  lines,  where  the  potential  values  to  be 
plotted  are  equally  spaced  values  of  V.  Also,  care  must  be  taken  when 
using  the  flux  line  plots  to  determine  the  regions  of  highest  electric 
field.  This  is  because  the  electric  flux  is  a volume  field  flow  which 
is  difficult  to  represent  clearly  in  a two-dimensional  plot.  Electric 
flux  distribution  is  shown  by  describing  many  tubes  of  flux,  each  of 
which  encloses  a given  amout  of  flux.  The  flux  density  is  greatest 
where  the  tubes  are  the  narrowest.  In  a two-dimensional  plot,  a tube 
may  appear  to  be  narrow,  when,  in  fact,  the  dimension  that  is  perpen- 
dicular to  the  plotting  surface  may  be  very  large.  Thus,  what  appears 
to  be  a high  flux  density  region  may  actually  be  a low  flux  density 
region.  Knowledge  of  the  three-dimensional  geometry  is  necessary  in 
order  to  interpret  a two-dimensional  plot. 

In  spite  of  the  difficulty  in  interpretation,  plots  of  electric 
flux  lines  are  sometimes  useful.  They  indicate  the  direction  of  the 
electric  field  more  clearly  than  equipotential  plots  and  can  also  be 
used  with  potential  plots  to  more  clearly  define  regions  of  high  elec- 
tric field  intensity  or  distortion. 

To  obtain  the  solution  to  Equation  35  for  the  purpose  of  plotting 
flux  lines,  the  computer  program  is  executed  exactly  as  for  the  poten- 
tial solution.  All  that  needs  to  be  changed  are  the  boundary  conditions 
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and  dielectric  constants  such  that  the  equal  U-value  lines  are  ortho- 
gonal to  the  equipotential  lines.  The  potential  function,  V(r,z),  has  a 
gradient  that  is  perpendicular  to  lines  of  equal  V at  each  point.  The 
negative  of  the  gradient  of  V is  the  electric  field  vector  function,  E, 
as  in  Equation  36. 

E = - VV  (36) 

For  the  function  U(r,z)  there  can  be  defined  a corresponding  vector 
function,  F,  given  by  Equation  37. 

F = - VU  (37) 

F is  perpendicular  to  the  lines  of  equal  U at  each  point.  The  requirement 
that  lines  of  equal  V be  perpendicular  to  lines  of  equal  U is  equivalent 
to  requiring  that  vector  function  E be  perpendicular  to  vector  function  F. 

The  requirement  of  having  E perpendicular  to  F can  be  analyzed  at  a 
boundary  between  two  dielectrics  as  shown  in  Figure  19.  In  (a)  the 
dielectric  constants  and  C2  are  normal  dielectric  constants  as  they 
would  be  used  in  the  solution  for  the  potential  distribution,  V.  In  (b) 
the  dielectric  constants,  e|  and  e^,  are  the  values  as  they  would  be 
used  to  obtain  the  U distribution,  n and  t are,  respectively,  unit 
normal  and  unit  tangent  vectors  to  the  boundary  between  the  two  media. 

Since  there  is  no  current  flowing  along  the  boundary  between  the 
dielectrics,  the  tangential  electric  fields  on  each  side  of  the  boundary 
must  be  equal.  This  condition  is  expressed  mathematically  in  Equa- 
tion 38. 

t.E^  = t»E^  (38) 

Also,  the  normal  component  of  the  displacement  flux  density  must  be 
continuous  across  the  dielectric  boundary  as  expressed  in  Equation  39. 

ei(n»Ei)  = E2(n»E2)  (39) 

Further,  F corresponds  to  E,  and  both  are  actually  solutions  to  the  same 
equations  with  only  the  boundary  conditions  differing.  Therefore,  F 
must  also  meet  conditions  that  correspond  to  Equations  38  and  39  for  E. 
These  conditions  for  F are  expressed  in  Equations  40  and  41. 


t.Fi 

= t.F- 

(40) 

G|(n*F^) 

= e^(n*F2) 

(41) 
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(a)  V problem. 


(b)  U problem. 


Figure  19.  Boundary  between  two  dielectric  media. 


E and  F can  be  expressed  in  terms  of  their  tangential  and  normal 
components  at  the  boundary  between  tlie  two  dielectrics.  In  medium  1 at 
the  boundary: 


*^1 

II 

rr> 

+ n.Ej 

(42) 

^1 

= f.F^ 

+ n»F^ 

(A3) 

In  medium  2 at  the  boundary: 

= t.E2 

+ n»E2 

(44) 

^’2 

= t.F2 

+ n*F2 

(45) 

The  requirement  that  E and  F be  perpendicular  must  hold  at 
boundary.  The  nonzero  vectors  E and  F are  perpendicular  if 
their  dot  product  is  zero  as  in  Equation  46. 

the  dielectric 
and  only  if 

E»F  = 

0 

(46) 

Substituting  Equations  42  and  43  into  Equation  46  gives  the 
for  perpendicularity  at  the  boundary  in  medium  1: 

condition 

Ej-Fj  - (£»E^ 

)(t*F^) 

+ (n»Ej^)  (n»Fj^)  = 0 

(47) 

The  condition  for  medium  2 is 
45  into  Equation  46  to  yield: 

obtained 

by  substituting  Equations  44  and 

^2*^2  ~ (t»E2 

)(t«F2) 

+ (n»E2)(n*F2)  = 0 

(48) 

Equations  47  and  48  are  equivalent  and  can  be  compared  by  substitution 
of  Equations  38,  39,  40,  and  41  into  Equation  47.  Equation  49  is  Equa- 
tion 47  written  in  terms  of  E2  and  F2 . 

G G * 

(t»E2)(t*F2)  + (n*E2)  ^ (n»F2)  = 0 (49) 

Equations  48  and  49  can  both  be  true  only  if  the  dielectric  constants 
satisfy  Equation  50. 
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(50) 


E 


£ 


2 

1 


Equation  50  holds  only  if,  for  any  value  of  a. 


e 


1 


(51) 


(52) 


For  simplicity,  a can  be  chosen  to  be  unity.  Equations  53  and  54  give 
the  dielectric  constant  changes  used  when  it  is  desired  to  obtain  the  U 
solution  for  plotting  electric  flux  lines  that  are  perpendicular  to 
electric  potential  lines  for  a given  dielectric  configuration. 


e 


1 


1 


(53) 


(54) 


In  addition  to  changing  the  dielectric  constants,  the  boundary 
conditions  must  also  be  chosen  so  as  to  give  a solution  of  electric  flux 
lines  that  are  perpendicular  to  equipotential  lines.  If  a grid  boundary 
corresponds  to  what  is  known  to  be  a flux  line,  that  portion  of  the  grid 
boundary  may  be  assigned  a constant  value,  U]^.  This  condition  often 
applies  to  the  section  of  the  axis  of  symmetry  between  the  electrodes  of 
an  insulator.  For  boundaries  that  are  distant  from  the  insulator  assembly, 
a curve  that  approximates  a known  electric  flux  line  can  be  estimated  as 
described  earlier  for  the  equipotential  solution.  The  points  along  this 
curve  may  be  assigned  a second  constant  value,  U2.  Some  portions  of  the 
grid  boundary  may  correspond  to  an  equipotential  line  or  an  electrode  or 
other  conductor.  The  electric  flux  lines  must  intersect  these  boundaries 
in  a perpendicular  manner.  This  condition  is  met  by  describing  the 
region  beyond  the  equipotential  line  or  electrode  boundary  as  having  a 
dielectric  constant  of  zero.  The  argument  here  is  consistent  with  that 
for  the  potential  distribution  problem  where  an  electric  flux  line  and 
zero  dielectric  constant  are  used  to  satisfy  certain  boundary  conditions. 
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Figure  20  shows  the  boundary  conditions  for  a potential  distribu- 
tion problem.  The  corresponding  boundary  conditions  for  obtaining  the 
electric  flux  line  solution  are  shown  in  Figure  21.  The  geometry  of  the 
problem  and  the  dielectrics  involved  are  designed  to  give  a simple 
illustration  of  the  concept  and  represent  no  actual  device. 

The  enlarged  equipotential  plot  is  given  in  Figure  22  while  both 
the  equipotential  lines  and  the  lines  of  displacement  flux  density  are 
shown  in  Figure  23.  They  must  be  referred  to  as  lines  of  displacement 
flux  density,  D,  since  their  density  is  continuous  across  the  dielectric 
boundaries.  It  is  more  practical  to  plot  the  flux  as  lines  of  D rather 
than  as  lines  of  E,  since  E is  not  continuous  across  the  boundaries 
between  dielectrics.  If  lines  of  E were  plotted,  they  would  need  to 
have  a density  in  the  dielectrics  that  is  inversely  proportional  to  each 
dielectric  constant.  Thus,  there  would  be  fewer  lines  in  the  higher 
dielectrics  than  are  shown  in  Figure  23. 

In  order  to  interpret  the  displacement  flux  density  plot,  it  is 
important  to  understand  the  meaning  of  the  line  spacing.  In  Figure  23 
the  line  spacing  is  determined  according  to  the  displacement  flux  density 
which  flows  in  the  plane  of  the  figure.  Therefore,  it  gives  only  a two- 
dimensional  representation  of  the  flux  distribution.  The  third  dimension 
must  be  visualized  by  rotation  of  Figure  23  about  the  axis  of  symmetry. 
The  significance  of  this  rotation  is  that  the  lines  nearer  the  axis 
sweep  out  a smaller  volume  than  do  the  lines  that  are  far  from  the  axis 
of  symmetry.  A given  amount  of  flux  in  a smaller  volume  has  a higher 
flux  density.  Therefore,  the  regions  near  the  axis  of  symmetry  have  a 
higher  relative  flux  density  than  appears  in  Figure  23. 


DISCUSSION 

The  computer  program  is  designed  to  yield  a plot  of  equipotential 
lines.  From  the  equipotential  plot  it  is  easy  to  visually  determine  the 
regions  of  highest  electric  field.  These  regions  are  characterized  by 
closely  packed  equipotential  lines.  In  the  equipotential  plot  there  is 
a fixed  potential  difference  between  any  two  adjacent  equipotential 
lines.  A relatively  short  distance  between  lines  Indicates  a high 
voltage  gradient.  As  discussed  earlier,  the  voltage  gradient  is  the 
dominant  component  of  the  electric  field  near  the  insulator  and  electrodes 
at  VLF  and  lower  frequencies.  The  regions  having  highest  voltage  gra- 
dients are,  therefore,  most  susceptible  to  electrical  breakdown. 

Since  the  equations  used  in  the  development  of  the  computer  program 
are  valid  for  all  frequencies  up  to  and  including  VLF,  the  program  has  a 
wide  range  of  application.  Most  notably,  it  is  a powerful  tool  that  can 
be  used  in  the  analysis  of  most  60-Hertz  power  system  insulators. 
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Figure  23.  Enlarged  view  of  equipotential  and  flux  line 
plot  for  problem  of  Figures  20  and  21. 
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Poisson's  Equation 


In  determining  ttie  potential  distribution,  the  program  essentially 
calculates  the  solution  to  Laplace's  Equation  20.  The  program  can 
easily  be  adapted  to  solve  Poisson's  Equation  55,  which  applies  to  prob- 
lems involving  a free  charge  distribution. 

= ■£  (55) 


In  this  case  the  charge  density,  p,  must  be  integrated  over  each  volume 
element,  and  Equation  29  becomes: 


vol 


dv 


^ D ‘ds 


s 


s 


vol 


(56) 


This  has  the  effect  of  replacing  the  zero  on  the  right-hand  side  of 
Equations  31  and  34  with  a term  that  represents  the  net  charge  distrib- 
uted throughout  the  volume  element.  Otherwise,  the  method  of  solution 
is  the  same  as  described  earlier. 

Other  Applications 

The  program  is  suitable  for  solving  many  other  types  of  problems 
that  involve  Laplace's  or  Poisson's  equation.  For  example.  Equation  57 
is  Laplace's  equation  as  it  applies  to  steady-state  conduction  heat 
transfer  in  regions  where  there  is  no  heat  generation. 

V^u  = 0 (57) 


Here,  u represents  the  heat  energy  or  temperature  as  a function  of 
position  in  space. 

In  problems  dealing  with  mechanical  stress  there  sometimes  arises 
an  equation  of  the  form  of  Equation  58. 

V^A  = 0 (58) 

This  can  be  solved  by  the  computer  program  through  a two-step  process  if 
there  can  be  found  an  X such  that 

V^A  = X (59) 
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If  the  boundary  conditions  for  X are  known,  the  complete  solution  for  X 
can  be  found  by  solving 


V^X  = 0 


(60) 


Equation  60,  which  is  Laplace's  equation,  is  identical  to  Equation  58 
with  the  substitution  for  A made  as  in  Equation  59.  Then  A itself  can 
be  found  by  solution  of  Equation  59,  which  is  Poisson's  equation,  uti- 
lizing the  boundary  conditions  for  A and  the  solution  for  X. 

Two-Dimensional  Geometries 

If  the  geometry  of  a problem  has  only  two  dimensions,  the  equations 
are  much  simpler  than  for  three  dimensions  with  one  axis  of  symmetry. 

The  volume  element  of  Figure  4 is  reduced  to  a surface  element  as  shown 
in  Figure  24  for  a rectangular  coordinate  system.  Instead  of  summing 
the  flux  crossing  surface  areas,  the  sum  of  the  flux  which  crosses  the 
boundaries  of  the  surface  in  Figure  24  is  taken.  Equation  61  is  the  new 
equation  which  corresponds  to  Equation  31. 


h.-, 


+ 


^8 


+ (^2  ^ "3  ^3^ 


1^  ^^4  h 


"5 


ir  ^"^6  S 

4 


"7 


— (e,  h^  + Eg  h^)  + — (e^  h,  + h^) 


h^  ''"1  2 


h^  ""2  “1 


'3  3' 


“fT  ^^4  ^2  ^5  ^^4^  h~  ^^6  ^3  ^1  ^1^ 

3 4 


= 0 


(61) 


Equation  61  is  a numerical  form  of  Equation  29  where  the  surfaces  have 
been  reduced  to  straight-line  boundaries,  and  a factor  of  one-half  has 
been  removed  from  each  term. 
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Appendix 


FVSOLVR  and  FVPLOT 


The  computer  listing  that  follows  is  given  in  the  form  of  two  main 
programs.  The  first,  FVSOLVR,  reads  input  data,  formulates,  and  solves 
the  equations,  giving  the  values  of  potentials  at  each  grid  intersection 
point.  The  second  program,  FVPLOT,  uses  the  potential  distribution 
solution  from  FVSOLVR  to  plot  lines  of  equal  potential.  It  is  convenient 
to  separate  the  process  into  two  programs  so  that  the  plotting  is  detached 
from  the  solution  calculation.  Thus,  it  is  easy  for  the  user  to  obtain 
several  different  plots  for  a given  problem  by  specifying  various  equl- 
potential  value  spacings  or  various  regions  of  the  grid  to  be  plotted. 
FVPLOT  can  plot  equal  value  lines  from  up  to  two  sets  of  data,  such  as 
equipotential  lines  and  flux  lines,  on  one  grid.  If  both  equipotential 
and  flux  lines  are  to  be  plotted,  then  two  solutions  from  FVSOLVR  are 
required  before  FVPLOT  is  used. 

FVPLOT  is  designed  to  produce  Calcomp  pen  plots  on  drum  plotters 
having  a width  of  11  Inches.  In  the  form  listed,  FVPLOT  utilizes  the 
CDC  7600  computer  at  the  Lawrence  Berkeley  Laboratory  via  remote  ter- 
minal access. 

Comment  statements  given  throughout  the  listings  indicate  the 
general  organization  of  the  programs.  A detailed  users  manual  is  in 
preparation. 


u>ro«*^o^  OD  ^<^u><p‘Uiru^o 


'! 


c 

c IKSOUATOR  field  valine  solver  (AXIAL  SYMMETRY) 

C 

DIMENSION  HX (100>  *hY (200) *HN(99) »HM(200>  iAN(99) t 
IEPSB(  199.  lo) .EPS^ (99. 10) .ANEPS (99) .RHN (99) .SYM (200*A) . 

22EPS1  (99)  .V  (100.2(^0)  'H  (200)  iAR  (200. 100) 

DIMENSION  Gh(200) .V0(200) 
integer  BCR(199.1w) •U02(99.10> 

integer  BLNN 

COMMON /aANAfiG/NDPN.NUM8LK.B.AR»NN.NN2.V.NX.NY.LX,LY.M.NCK8 
COMMON/HORARG/BOR.EPSP 

common/ verarg/bdz.epsz 

COMmOn/SPArg/GH, VG.HM 
DATA  oLNK/Am  / 

DATA  aR/?000O*0.O/.SsYM/8oo*0i0/.B/2O0*0.O/. V/2oOOO*"1 .0/ 

C 

c fohmai  statements 
C 

format  (loEa.o) 
format  (2oU) 

FORMA) ( 16F5.0) 

FOHMAI (8F10. 3) 

format (1H1.4HNX  s.Ia.AX.AHNY  a,l4) 

FCRMAKIH  .A1.3HHa(.I3.3H)  « .F8.3  . A ( AA  .6X  . 3HHX  ( . 1 3 . 3H)  3.F8.3)) 
FOHMATdH  .A1.3HHY(»I3t3H)  « »F  8 . 3 .4  ( A*  .6X  . 3HH  Y ( . 1 3 » 3H ) a.Fe.3)) 
FORMAT (1H1.9X.1H..14F8.3) 
format (loX.lH*.14Fe«3) 

FORMAT  ( 10*»1H*.  14  (8M—-*  — ) ) 
format (1X.F8.3.2H  ♦.14F8.3) 

FORMATTlHO.lSHlNKIAL  vALuE  ■.Fe.3) 

FORMAT (1H0.6HNPTS  *.I7) 

FORMA)  (1X.7)-NUMBLN«.15) 

format  (IX. 5t-NCK5«. 15) 
format (1X.6HNTEST*. IS) 
format (13H  nONSYMmETRIC) 
format  (IHI) 
format (IH  ) 

FORMA) (IHO) 

FO^HAI (IH*) 

format  (iH  .27HREADY  TO  BEGIN  CALCULATIONS) 

FORMA) (2I8.4F12.3) 

C 

C initialization 

c 

read  2 .NX.NY.NEPRtNEPZ.IFBOY 

print  5 .NX. NY 

print  2o 

NXMISnX-1 

nymi*ny-i 
nn*nx 
NN2*Nn*2 
NNPl»NN»l 

rewind  lo 
PEwlNO  33 
rewind  34 

c 

c read  grid  data 


BES)  AVAIUBLE  CG?i 


47 


o rt  o o 


c 


r 


REAU  1 «(HX(I)«I>ltNX) 

KlsQ 

NCOL»N*M1/5o 
eCTO  120 
loo  NC0L»NC0L-1 
UO  INK»NCOL*50 
130  M«K1*1 
^2*K1 ♦ INK 

PRINT  6 f (BLNKtK.H*(K) tK=KlfK2ibO) 

IE  (K1.EO.SO>OR.K1.Eu.NX)  goto  IAO 
IF  (Kij.tQ.N*)  100*130 
lAO  read  1 * (my (J) * J*1*NY) 

PRINT  18 
K1  = 0 

NCOL*NYMi/5o 
GOTO  l70 
150  NCOL»NCOL-l 
170  INK«NcOU*50 

IHO  M = Kl*l 

K2sKl*INK 

PRINT  7 * (eLNK*KtHY (K> iK=K1*K2*50) 

IF  (Kl.Ea.bO.OR.Kl.Eu.NY)  goto  200 
IF  <K2.E(3.Nt)  150*1S0 

2t'0  read  2 t ( (dOR  (I*  J)  *J*1  tNEPR)  *I«l*NYMl» 
220  read  2 tKO,KN.KK 

READ  3 t (EPSR (KO*K) (KslfKK) 

IF  (K0.E0*KN)  goto  260 
KISKO+I 

00  2A„  I»KI,KN 
DO  2Ao  K»1*KK 
240  EPSR(IiK)aEpSR(KO»K) 

260  IF  (Kk.EQ.NePR)  GOTO  300 
^PSRK♦1 

DO  28o  I*KO»KN 
DO  28j  J»KP,NEPR 
280  EPSR(I*J)«EpSR(I*KK) 

300  IF  (KN.LT.NYM1)  GOTO  220 

read  2 *((BOZ(I*J)*J»l»NEP2)*Ial»NXMl) 
320  read  2 tKO,KN*KK 

read  3 * (£PSZ(kO*K) *K«l*KK) 

IF  (Ku.EO.RN)  GOTO  360 
M*KO*l 

DO  34u  I«KI,RN 
00  34j  K»1*rK 
340  EPSZ(I*K)*EPSZ(K0.K) 

360  IF  (Kr.EO.NePZ)  goto  400 
KP*KK*1 

DO  38u  IsKO.KN 
DO  3Bo  JaKP.NEPZ 
380  EPSZdf J)>EpSZ(I*KK) 

400  IF  (KN.LT.NXM1)  GOTO  320 
IF  (IFBDY.EO.O)  GOTO  550 

INTERPOLATE  fine  GRID  BOUNDARY  POTENTIALS 
FROM  COARSE  GRID  POTENTIAL  SOLUTION 
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PEAO  2 iNGX.NGY 

REAO  4 t (GH(I) tlaliNGX) 

read  4 f (Ve(I) tI»l»NGX) 

CAUL  SPCOEF  (NGX) 

DO  450  UltNX 

450  V(I.Nt)»SPLINE(NGXtHA(I) > 
read  4 . (GH(I) fl-ltNGV) 
read  4 t (V6(I) tl"! ‘NGY) 

CALL  SPCOEF  (NGY) 

DO  500  I»1*nYm1 
5u0  V(NX.J)»SPLINE(NGY*HY(1)  ) 

C 

C read  dOONOARY  CONOITIONS 

c 

550  read  2 fNROTS 
PRINT  le 
00  GOu  N«ltNPOTS 
read  1 .POT 
print  12.PCT 

read  2 .nlines 

PRINT  2 .NlINES 
00  600  K»1»NLINES 
read  2 »L1,L2.m1.M2 
PRINT  2 .Ll,L2»MltM2 
DO  60w  I»LI,L2 
DO  600  J«Ki,M2 
600  Vd.JisPOT 
C 

c Parameter  imtialuation 
c 

NUMBLk^O 

NCK»0 

NCKl»l 

NCKa*NN 

NCK5»u 

NCNT» j 

LY»2 

Z=«MY(i)*hY(2) )/2. 

2PHa(MY(2)  .(hYO)  )/2. 

MjMX*NX-1 

IF  (MX ( 1 ) .NE.O.O)  NUM*aNX-2 
NPTS»NUMX* (nY-2) 
print  20 

PRINT  13.NPTS 
nCK7®nPTS-NlMX 
NO  I AGs  I 
NCPl*NOlAG.l 
NOPNsNDI AG*MJMX 
print  IB 
DO  65o  6«1.nXM1 
650  2EP51 (K)«EPSZ(K.1) 

C 

C CALCOlAIE  distances  bET.<EEN  GRID  LINES 
C 

DO  700  J*1.NXM1 
7o0  ( J) =HX { j* j ) -HX < J) 

DC  72i,  Jsl,^yMl 
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BEST  AVAIUBlt  COPY 


r 


7<.0  Hf'(J)sHTr  (^♦l)-HY(j) 

c 

C CftLC^JLAlE  AkEaS  of  hurIZOMAU  faces,  set  rhn  array 

c 

F»HX ( I ) 

AN  (1)S0,5*MN (1 )*(K*HN(1)/4,) 

RHN  ( 1 ) sHthN ( 1 ) /2. 

00  75a  J»2,nXM1 
HhN(J)*(hX(j)*HX(j*1>)/2. 

750  AN  (J)a(.HHN(j)**2-riHNi(J-l)**2)/2. 

C 

c bCTTON.  Eoee  area  - ers  oata  set-up 
c 

EPSl=tPSR(l,l) 

R = hX  < l) 

LX=1 

J=1 

CALL  fiOHiZ  ( 1 , JfEPSl  ,H,HrtN  (LX)  ,ANEPS  (LX)  , AN  (LX)  ,LX,HX  (LX)  ) 

DO  aOu  LXa2,N|XMl 
J=1 

CALL  hOKlZ ( 1 ,J,EPS1 »RHN(LX-1) , RHN (LX) » ANEPS (lX ) , aN ( LX ) ,lX ,hX ( LX ) ) 
eyo  CCNTInUE 
PRINT  22 

print  20 
(iOTO  3000 
C 

C SHIFT  BLOCK  UP  ANU  ZERO 

C 

'SOO  00  1000  K»1,NN 
KPN^a^♦^^ 

DO  95a  I»l»4 

SYM(K,I)«SYh(KPNN,I) 

950  SYM<KPNN,I)aO.O 
B (KPNN)*0,0 
DC  lOi/O  1*1, NN 
1000  AR (KPnN.I ) »0t0 
C 

C CALCULATE  BLOCK  OF  aH,  B ARRAY  VALUES 

C 

3000  DO  36i)0  L*NnP1,NN2 

IF  (NCK,E(J.q)  GOTO  3010 
lx*lx»i 

I R=HX(lX) 

I IF  (Lx.lT.NX)  goto  3400 

\ LY=L7*1 

1 IF  (LY.EQ.NY)  goto  4000 

1 Z*ZPH 

: ZPHsMY (LY) *hM(LY)/2. 

3010  J=1 

NCKaNCK*! 

EPSiatPSR(LY,l) 

IF  (MXd)  .NE,0*0)  goto  3200 
lX*! 

R*0.0 

IF  (V (LX,LY) .GE.OtO)  GOTO  3300 
i SYM(L,1)*0.0 

I SYM(L,3)aANEPS(l)/HM(LY-l) 

BEST'AVAIIABIE  COPY 

1 


L 


I 


1 

I 


I 
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CALL  mOHIZ (lYi J*Et^Sl |R»«MN (L*) iANEPS(L*)  »aN(L*)*LX*mX(lH)) 
SYM(L.4)»ANEPS(1)/HM(LY) 

CALL  vEhT (LXf ZEPSl (LX) , Z i ZPH t RH\ ( L X ) .AZEPStLYtHY(LY)) 
SYM(L.2)*AZePS/HN(LX) 

AR(L*iNDlAG)s-SYM(Hl)-SYM(Lt2)  -SYM(L«3)-SYM«Lf4) 

IF  (V(LAiLY-l) .LT.O.Q)  goto  3170 
a(L)=-V(LXtLY-l)*bYM(i.t3) 

SYM(Lt3)»0.0 

3170  IF  (V  (UX*1  ,t.Y)  .LT  .0*0)  GOTO  3180 
B(L)=3(L)-V(LX*ltLYr*SYM(L.2> 

SYM(Lf2)*0.0 

3100  IF  (V(lX,ly*1)  .LT.OjO)  goto  31'J0 
b(L)=el(L)-V(LX*LY  + l)*SYM(L»<>) 

SY'l(Lt4)a0.0 

3190  AR  (HinOP)  ) =SYm  (Lt2) 

AR  (Lti'iOPN)  =SYM  (L*<f) 

GOTO  3600 
3200  LX=2 

R = mX  (lX) 

IF  (V (LX.LY) .GE.O.O)  GOTO  3330 
SYM(Li 3) =ANtPS (2) /HM  (LY-1 ) 

CALL  nORl2(LY»JtEPSl tMHN (LX-l ) tRMN (LX),ANEPS(LX)*AN(LX)tLX,HX(LX)) 
SYM(l.4)=An6PS(2)/hM(lY) 

CALL  yEWT(lx-H2EPS1  (LX-l  ) t Z t ZPh  t WHN  ( LX- 1 ) f AZEPS  . L Y i HY  ( L Y ) ) 

SYM(Lil )=AZEPS/hN(LX-1) 

call  vEHT  (LX.ZEPSI  (LX)  tZ.ZPHiHhiN  (LX)  . AZEPS  *LYtHY  (LY)  ) 
SYM(L.2>=AZEPS/HN(LX) 

Afi(L»i'iOlAG)s-SYM{L»l  >-SYM(L*2)-SYm{L*3)"SYM(l.4) 
fc((L)=-Y(l tLY)*SYM(L*l) 

SYM(L.1)=0.0 

IF  (V  (LX,ly-1 ) .LT.y.O)  GOTO  3270 
b(L)=d(L)-V(LXtLY-l)*SYMlL*3) 

SYM(Li3)=0.0 

3270  IF  (Y(LX*1,lY)  .LT.^j.O)  GOTO  3280 
b(L)»8(L)-V(LX»ltLY)*SYM(t»2) 

SYi«1(L.2)=0.0 

32«0  IF  (Y  (LX,LY*1 ) .LT.0?0)  GOTO  3290 

b(L)=d(L)-V(LX«LY*l)*SYM(L*A> 

SYM (L i4) »0 .0 

3290  AR(Lt\0Pl)=SYM(L*2) 

AR  (lh'.OPN)  aSYP  (L*A) 

GOTO  3600 

3300  AP(L*‘ND1AG)=1. 
b(L) =V (LXf lY) 

call  riOHIZ(LY.J»EPSltR,RHN(LX)  ,Al\(EPS(LX)»AN(LX),LXfHX(LX)) 

CALL  vEHT  (LXfZEPSl  (LX)  t Z . ZPH  i RHN  ( L X ) ,AZEPS*I-Y,hY(LY)) 

GCTO  36o0 

3330  Afi  (UNOI  AG)  =1  , 

B (L) =v (LXtLY) 

CALL  30HlZ(LY,JtEPSl.PHN(LX-l) ,RMN(LX) ,AN£PS(LX) ,AN(LX)iLX,HX(lX)) 
CALL  yERT(lXiZEPSi (LX) t Z » ZPM » RHN ( L X ) » A ZEPS * L Y » H Y ( L Y ) ) 

GOTO  3600 

3400  IF  (V(LXtLY) .GE.O.u)  GOTO  3330 
SYM(L.3)=ANEPS(LX)/HM(LY-1) 

CALL  r.Or(IZ(LYtJfEPSl,RHN(LX-l)  »RHNilx)  ,ANEPb(Lx)iAN(Lx)»Lx,HX(LX)) 
5YM(Lt4) sANEPS (LX) /HM(LY) 

SYM(L.l )»AZePS/mN(LX-1 ) 
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BEST  AVAIUBIE  COPY 


rj  r)  o r>  r>  o 


CALL  vEt'T  (UXtiEPSl  (LX)  tZtZPHtRMN  (LX)  . AZEPs  • L Y • h Y ( L Y ) ) 
SYM{L.2)»AZtPS/HN(L)^) 

AC(L*^DIAe)a-SYM(t»l)-SYM(l.l2)-SYM(L*3)-SYM(L»A) 

IF  (V (LX-l ,lY) .LT.0.0)  GOTO  3460 
B(L)=-V(L*-liLY)*SYM(L»l) 

SYM(l.1)«U.O 

3460  IF  (V (LX.LY-l) .LT.0*0)  GOTO  3470 
B(L)*b (L)-V {LX*LY-1 )»SYM(Lt3> 

SYM(L*3)»0.0 

347o  IF  <v (LX*1*ly) •LT.u.u)  GOTO  3460 
b(L)»6(L)-v(LX*ltLYy*SYM(Lt2) 

SYm(Li2>*0.0 

3460  IF  (V (LX,Ly»1) .LT.0.0)  GOjO  3490 
B(L)*d(L)-vlLX»LY4iy*SYM(Lt4> 

SYM(Lt4)a0.C 

3490  AR(Lfl'.0Pl)sSYM(Lt2) 

AR  (HM)PN)  «SYH  (L  !<♦) 

3600  continue 
4000  NTEST=3 
C 

C »«fiITE  BLOCK  OF  AR|  b ARRAY  VALUES  ONTO  TaRE 

C 

WRITE  (34)  (b(N) * (Art(NtR) ,MaltN0PN) *NaNNPl.NN2) 

NTe.STa4 

NUMHLKaNUWSuK* 1 
NCNTaNUMeLK*NN 
IF  (NUniSLK.EC<1)GOIO  4300 

CHECK  MATRIX  symmetry 

4100  00  42j0  L»NCK1»NCK4 
NCK5=nCK5*1 
NCK2»l*1 

NCK3at*NUHX 

IF  (SYM(Lt2) ,NE.SYM(ncK2* 1 > y GOTO  4150 

IF (nCk5*GT .NCK7)GOTO  4200 
IF  (SymIl»4)  .E(J»SYM(NCK3»3)  ) GOTO  4200 
4150  print  17 

print  23»NLM0LK,NCK5,SYM(NCK2tl) ,SYM<l,2) #SYM(NCK3»3) »SYM(L*4) 
4200  continue 

IF (NCKl .EC.NNPI )GOTO  4500 
4300  IF (NCnT;lT.nPTS)GOTO  900 
NCKlaNNPl 

NCK4*NNP1 ♦NPTS-NCK5-1 
GOTO  4100 
4500  continue 

PRINT  15«NCK5 
NCKHaNPTS-(NUMBLK-l)*NN 
print  2 tNCKS 

ZERO  aR  and  8 arrays 

DO  4800  J>ltNN2 
B( J)*u.o 
DO  48C0  L*ltNN 
4800  AR(JfU)*0>0 
C 
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AVAIIABLE  copy 


ooooo  o no  o r>  o o oo 


C 

C 


CALL  dANOEl)  MaTHIX  solver 

IF  (MXH)  .NE»0»0)  M»2 
CALL  dANSOL 
NTEST«5 

print  UtNTEST 
print  2 tLX.LV 

WRITE  SOLOTION  ONTO  TAPE 

wPlTE(10t2)  NX, NY 

WRITE (I  0,4)  (HX(J) tJxliNX) , (HY(K> ,K*1,NY) 
WRITE  (I 4)  ( (V (I, J) ,I>1,NX) • JsIfNY) 

END  FILEIO 

print  SOLLTION 


5000  M2so 


Ml»5o 

5100  M2sm2*5o 

IF  (M2.LE.NY) 
M1«NY*50-M2 


GOTO  5200 


M2»NY 
N2»0 

5300  M*N2*l 
N2*N2*14 

IF  (N2.GT.NX)  N2»NX 
PRINT  a , (hiX(K)  ,N»N1,N2) 

print  10 

DO  5<H)0  J»l,Ml 
IsM2*i-J 

00  PRINT  V1»MY (I) * (V(N,I ) ,K»N1 »N2) 

print  10 

print  9 , (hX(K) ,n»N1,N2) 

IF  (N2.LT. NX)  GOTO  5300 
IF  (M2.LT.NY)  GOTO  5100 
STOP 
END 

SUHROgTiNE  dANSOL 


SOLVES  8AN0E0,  SYMMETRIC  MATRIX 

common  /BANaRG/  MM,NUM8LK,B»A,NN,NN2, V,NX,NV,LX,LY,LM8,NCKp 
DIMENSION  d(200)*A<<JOO*100)*V(100,200) 

NL»NN. 1 

NHxNN.NN 

LXxnX 

LY»NY-1 

REwINu  34 

N6*0 

GO  TO  150 


REOUCfc  E(3UAtI0NS  dY  dLOCNS 


I.  SMiFf  dtOCN  OF  EOUaTIONS 
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or>r>  ooo  ooo  r)or> 


C 


100 

DC  1?3 

^^SMN«N 
tt  {N»  «tj  (NM) 

B(\'M)so.n 
DO  12d 
A (NtM)  aA 

125  A(NMtM)so.O 

2.  flEriD  next  block  0^  EUUATIOnS  INTO  COHE 

IF  (NyMBUK-NB)  l5j»200*l5j 
150  PEAD  13<*)  (B  (N)  » <A  (iNty)  *M*1  tMH)  tN*NL  •NH) 

IF  (Nd)  200»l00*2o0 

3.  bEljdce  block  of  equations 

200  DO  300  N*lfNN 

IF  (A<N»i))  225*300*225 
225  B(N)=d(N)/A(Ntl> 

DC  2Td  L»2*C'M 
IF  (A(N.l))  230*275*230 
230  C*A(N,L)/A(^,l> 

I«N*L-1 
J = 0 

do  25u  K*L|^M 
J»J*1 

25o  A(I*J)»A(I,J)-C*A(N*^) 

d(n»d<I»'AtN*L>«d(N) 

A (N,U) aC 
275  continue 
300  continue 

4.  a«iTE  Block  of  heouceo  equations  on  Tape  2 

IF  (NoMriLK-NB)  375*A00*375 
375  WRITE  (33)  IB  (N)  * (A  (N*0)  *Ma2tM(*i)  *Nal,NN) 

GO  TO  luo 


back-substitution 


4o0  00  45y  '^al*NN 
Ti*NN*  1-M 
DO  A2d  Ka2,^M 
LaN*K-l 

425  B<N)«ri(N)-A(N|K)*olL) 

N^<BN♦^N 
B (NM) aB (N) 

IF  (Nri.EO.NLMBLK.AND.N.GT.NCKB)  GOTO  450 
LX»L*-1 

IF  (LX.GE*LBB)  goto  445 

LY-LT-1 

LXbnX-I 

445  V(LX*LY)aBIM 
450  CONTINUE 
NB*NB-1 


BEST  AVAIWBIE  (DPT 


o or»  rvo  r»o  rir»r>r> 


IF  (Nd)  475t500»*75 
475  backspace  33 

read  <33)  (8  *Ms2tMP)  ,N»i  ,NM 

BACKSPACE  33 
GO  TO  400 
500  CCNTLmOE 

return 

END 

subroutine  HORIZ<UYf J,EPSi tRA,RBtANERStAREA*LX,«) 

calculates  dieleciric-area  product 
OF  HORIZONTAL  FACE  OF  VOLUME  ELEMENT 

dimension  EPS(199«10) 
integer  80R(l99tlu) 

COMMON/HORARG/BDRtEPS 
IF  (BuR(LT»u>-2*LA)  l,2t3 

1 J»J*1 

EPSl«tPs<LY,J) 

IF  (BU«<LTtj»-2«LA)  l,2»3 

2 J»J*1 

EPS2»tPS(LY,J) 

ANEPS»0*5« ( (EPS1-EPS2)»R**2*EPS2*RH**2-EPS1*RA**2) 

EPSl«cPS2 

return 

3 ANEP5»ePSl*AREA 

return 

END 

subroutine  vert (LAtEPSl *Z.ZPH*RfAZEPS.LY*Y) 

CALCUlAIES  uIELECIRIC-AREA  product 

OF  vertical  face  of  volume  element 

dimension  EPS(99*10> 
iNTEGtR  BDZ(99»10) 
common/ VERaRG/BOZ I EPS 
K«l 

IF  (BuZ<L*fK)-2*LY)  1,2*3 

1 K*K*l 

EPSl*tPi»(LX,K) 

IF  (BcZiLX,K)-2*LY)  1*2*3 

2 K»K*1 

EPS2«tPS(LX,K) 

AZEPS»R*(EPSl*(Y-Z)*EPS2*(ZPM-Y)  ) 

EPSl«tPi>2 

return 

3 AZEPSaEPSl*R*<ZPH-Z) 

return 

END 

subroutine  SPCOEF  (h) 

calculates  COEFICIENT?  OF  CUBIC  SPLINE  USED  IN  INTERPOLATION 

DIMENSION  Xn<200) *EN(2oO) , S (2oo ) *RhO (2oO ) * TAU (200) 
COMMON/SPARe/XN*FN*S 


BEST  AVAIWBLE  COPY 

A 


no 


Nf'aaN-Z 
PHO(2)>j.O 
T«U(2)>U.O 
DO  I la2,NMl 

IPl-l*! 

HIM1«XN(1)-XN(IMI) 

MI«XN(IPl)-XN(I) 

Temp* (HiMi/hi ) • (Rho (I ) ♦2*) ♦Si 
PH0(I*1)«-1,/TEMP 

0»6.*( <FN(1P1)-FNII) )/HI-(FN(I)-FN(IMi) 

1 TflU<Itl)»(0-HIMl*(AD(I>/HI)/TEMP 
S(1)»0. 

S(N>«o. 

DO  2 lalfNM^ 

IB«N-I 

2 S(IB)«RH0(IB*U*S(IB^1)  tTAUdB^l) 
return 

END 

FUNCTION  spline  (N,X) 

evaluates  Spline  function 

DIMENSION  XN (200) tFN(200) tS(SOO) 

C0MM0N/SPAB6/XN,FNtS 
IF  (X.GE.Xn(I))  goto  1 
Hl«XN(2)-XN(l) 

SPLINe»FN(i)  ♦(X-Xn(1)  )*(  (FN(2)-FN(i)  )/w4-I-1*S(2)/6.) 

RETURN 

1 IF  (x,le.xn(n>)  goto  2 

NPl«N-l 

HNMI«XN(N)»XN(NMI) 

SPLINEaFN(N)^(X-XN(N) >•( (FN(N)-FN(NMI) )/HNMI4MNM1«S(NM1)/6.> 

return 

2 DO  3 I>2fN 

IF  (X.LE.XN(I))  goto  4 

3 continue 

4 L»I-1 
LPI»L4l 
A»XN(LPn-X 
BbX»Xn (L) 

HL*XN(LPl)»XN(L) 

SPLINE«A»S(l)«<A»«2/ML-ML)/6.^B«S(L*1)*(B*«2/HL-HL>/6, 

/ ♦(A*FN(L)4B*FN(LPI) )/ML 

return 

END 
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BESTAVAIUBIE  COPY 


FVPLOT 


C 

C contour  line  plotieh 

c 

OIMEN<jIUN  V (100«200>  (HX  (100>  *HY  (200)  ilBOF  (60 ) iGX  ( lOO  ) tOY  (ZoO  > i 
/OdOOiZOO)  tLXdSO)  iLYdSO) 

C0MM0N/0RAmaR6/X(600) fY(6U0>  *XMAXtYMAX •NQUAOtFiRSTX *DELTAX. FIRST Yt 
/DELTAY 

1 FORMAl (0FIO.3> 

2 format  (2014) 

3 FORMAT (6l4fF5.0) 

4 FORMAT (I4*FI0*0) 

5 format  (5M1NX  -flAtbXfAHNY  -tlA) 

6 FORMAT  (lhu»113tli4f I2Af II4) 

7 format  (7HoX  FR0M,F10.3*4H  TO»Fl0.3*eXt6HY  F ROM*F 10 *3 t 4H  TOt 
1F10,3) 

C 

C READ  potential  GRID  OATA 

C 

READ  4 flVUfUPOT 

read  3 (lOf INtjOtJNtNUUAOfKK.SZ 

read  3*  NX»NY 

READ  1»  (HX(I),  I*1iNX),  (HY(J)*  J»HNY) 

READ  It  ((V(ItJ)t  I«ltNX)t  J«ltNY) 

ISCANal 

IF  (IvU.EG.D  GOTO  ISO 
C 

C read  flux  GRID  data 

C 

READ  3 tRHiNHtKRA 
read  3t  MXiMY 

read  It  (GX(I)t  IsltMX)t  (6Y(J)t  JmltMY) 

READ  It  ((U(ItJ)t  I»ltMX)t  j«ltMY) 

C 

C shift  OATA  FOP  region  TO  BE  PLOTTED 
C 

IF  COtEOtl)  GOTO  120 
DO  lOo  I>IOiIN 
GX(I*1-I0)*GX(I) 

DO  lOO  J»ltHY 
100  U(I«l-IOtJ)sUdtJ) 

120  MX«IN*l-10 

IF  (J0tE(3tl)  GOTO  150 
DO  140  J«JOtJN 
GY(J*l-JO)«eY(J) 

DO  14o  iBltpX 
140  U(ltJ*l-jO)«U(ItJ) 

150  IF  (lUtEOtl)  GOTO  170 
DO  160  l>10tIN 
HX(I*l-IO)«HX(I) 

DO  160  U>ltNY 
160  V(I*l-10t0)»V(ltJ) 

170  NX»INtl-10 

IF  (JO.EOtl)  GOTO  200 
DO  100  JBUOtJN 
HY ( J*l-dO) »hY ( J) 

DC  ISu  iBltNX 

leo  v(itj^i-j0)=v(itj) 
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BEST  AVAILABLE  COPY 


a rtn  or»r»  o on 


2U0  NY«JN»1-J0 

PRINT  b *N;t*KY 

PPInT  t tlCtlNtJOfON 

PRINT  7 tHX(l) fHX(NX) ,MY(1) *HY(NY) 

X|MAXsHX(NX).HX(n 

YKAX«hV<NY)-HY(l) 

KS»0 

NXM1»nx-1 

NYM1«NY-1 

initialize  plot  and  draw  axes 

CALL  PLOTS  (Ht0»99) 

CALL  PACT0R(0«2S) 

CALL  PLOT  (o.»-U.»-3) 

CALL  HLOT(o.tO>5i-3) 

FlRSTXaMXIl) 

FIHSTy»MYI1) 

IF  (XmAX.GT.YMAX)  goto  220 

IF  (NLUA0.£(;.2*AND.  <2.*Xmax>  .GT.Ymax)  G0T0  220 
DELTAX«XmAX/1o. 

0ELTAY»DELTaX 
GCTO  24u 

220  OELTAY«YmAX/10* 

DELTAasDELTaY 
IF  (NgUADtEG.E)  GoTu  290 
240  IF  (NwoAO.Gl.l)  GOTO  2R0 
IF  (XmaX.GT.YPAX)  GOTO  260 
YAX»(YMAX/XwAX)*Io. 

CALL  AXIS(0.*10»*«:t«Z-AXIS  AXIS  OF  SYMMETRY, 4271 YaX,o.,FIRSTY, 

iDELTAYtu) 

CALL  AXlS(0.*10.tt)MW-AXISt-6tlO»»270.tFIRSTX,OELTAXi0> 

GOTO  290 

2b0  XAXa(AMAX/YwAX)*lu< 

call  AXIS(0.f0»*27HZ-AXIS  AXIS  OF  SYMpE Tr Y , ♦?? * 1 0 . »90 • *F IRS T Y , 
IDELTAY.o) 

CALL  AXIS(0^•0•  t6MWAXlS,-<,»XAXfO*»FIRSTX,DELTAXto) 

GOTO  290 

2e0  DELTAAa2.#i;tLTAX 
DELTAYa2,*OfcLTAY 
290  IF  <IvU.EG.2)  goto  305 

DRAW  tOOlPOTENTIAL  LINES 

DO  30u  N«ItNK 
CZ»K-nS 

z»cz*sz 

300  CALL  DRAW  (2tNXtNYfNXyiiNYMIthXtHY«ViISCAN> 

IF  (IvD.EG.J)  305f500 

calculate  and  plot  Flux  line  oisiribution 


r 


1 


pe*v(i,i) 

IF  (UP0T,GEtPA*AND*OP0T.LF.PU)  GOTO  320 
PA»Pb 
GOTO  310 
320  IP«I-1 

P»HY ( IM) ♦ (UPOT-PA)*(HY ( I)-HY (IM) ) / (Pfl-PA) 

H»MY (NH) -hY (MH) 

ISCAN32 

C«0. 

PAshX (KHA) 

330  ^P8■^HA♦1 

IF  (KkB.GT.nX)  GOTO  500 
RB*HX (K«B) 

EAa(V(KHAtNH)-V(KKAtMH) ) /H 
EB» ( V (KPBiNh) -V (KHBiMH) ) /H 
A2«(  (t8-EA)/  (pa-fiA) >/2. 

Ai3EA.PA«2.«A2 

CINC*u2* (Ra«*2-PA*»2) ♦Al*(R8-PA) 

IF  (R.LT.Pa)  CINCaA2*(«**2-RA**2)*Al*(R-RA) 
CsC*ClNC 

IF  (R.LE.Pd)  GOTO  3A0 

KPAsKhB 

PA>Ra 

GOTO  330 

3<>0  Z«U(XMAfMh)  ♦ (R-RA)  * (U(KRbttMH)  - U(^RA»MM)  ) / (RB-Ra) 
CALL  ORAW  (ZtKX,NYiNXMlfNYMl»GXtGY.U*ISCAN) 

3b0  RX*R 
CSOMs^. 

360  A2«( <tB-EA>/ (Rb-RA) )/2. 

A1«EA-RA*2.*a2 

CINCSU2* (Rb**2-RX**2) ♦A1*(rB-Rx) 

IF  (C-CSUP-CINC)  38v)t39ot370 
370  CSUMacSUM*CINC 
KRA3^HB 
KfiSaKHA+l 

IF  (KHB.GT.^X)  goto  500 
Ra»mb 

RX»RA 
PB*HX  (Kr(B> 

EA»(V(KRAtNH)-V(KHAtMH) >/H 
E0*(V (KRB»NP)-V (KRbtMH) )/M 

GCTO  360 

380  Ao»CSuM- (C*a2*R***2*AI*«*) 

CALL  KSOLY(N*flXfHB*A2*Al*A0) 

GOTO  400 
390  P»R8 

AGO  i;  = 0(KHA»MP>*  (R-RA>*  (0(KRM*MM)  -0(KRA»MH)  ) / (R8-Ra) 
CALL  ORAW  (2fNXfNYiNXBlfNYMlfGX«GY,U* ISCAN) 

IF  (IbCAN.EG.O)  goto  500 
IF  TR.NE.RtJ)  GOTO  350 
^PA»^r^B 
KPbsKnB^l 

IF  (KrtBtGTtpX)  GOTO  500 
Re»HX (KRB) 

PA«R 

GCTO  350 
C 


BEST  AVAIUBIE  CCri 
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c 

c 


OKAki  cLEcTHCOE  and  dielectric  CONUGURaTIcN 


1 

i 

t 


1 


t 


Soo  >^EAO  3 fNtlNES 
DO  70U  L»lfNLINES 
READ'  3 tNRTS 

head  2 ,<LX(N)tLY(N)iNal,NPTS) 

IF  (X^aX.GT, YMAX)  SOTO  60u 

IF  tNyUAO.E(;.2»AND.  (2.*XMAX)  .GT.YMAX)  gOTQ  600 

DO  5Au  Nal,NPTS 

LXNnN«lX(N)*i-IO 

LYNNNaLT  <N) ♦1-JO 

Y(N)xmY(LYNNN) 

GOTO  (Slct530«d30«520) •NOdAD 
blO  X (N)axMAX-r(X  (LXNNN) 

GOTO  S9%3 

520  Y (N) xYMAX-nr (LTNNN) 

5i0  X (NJxxMAXtHX (LXNNN) 

5aO  CCNTInOE 

X(NP7a+n»FiRSTX 
X «NPt5*2)»DELTAX 
Y(Nptb*l)«FIRSTY 
Y (NPTis  + 2)»DtLTAY 
CALL  line  (Y*X*NPISf lt0*0> 

GOTO  (7J0«5ti0t550f 550) iNUuAD 
550  DO  S6o  N«^,^PTS 
560  X {N)  »«!.*XPAX-X  <N) 

CALL  line  (YfX*NPTS»l»OtO) 

DO  S7j  N«1,nPTS 
570  Y (N)«i,*YPAX-Y (N) 

CALL  line  (YtXfNPIS*! fOtO) 

560  DC  59;j  N«l,^PTS 
590  X <N) *2»*XPAX-x (N) 

CALL  line  (Y*XtNPIStl*OtO) 

GOTO  700 

600  DO  6*u  N»l,^PTS 
LXNNNaL* (N) ♦l-IO 
LYNNNxlY (N) ♦1-JO 
Y(N)*riY(UYNNN) 

GOTO  (61o«630t520»520) tNODAD 
610  X (N)*MX(LXN^N) 

GCTO  6A0 

620  Y (N)*YHAX-hy (LTNNN) 

630  X (N)»XMAX^HX (lANNN) 

640  CONTlhOE 

X(NPT6^1)«FIRSTX 
X(NPt5*2)«0ELTAX 
Y(NpT5*1)»FIRSTY 
Y(NPts^2)»DELTAY 
CALLlI'^E  iXtYtNPIStltOfO) 

GOTO  (700«680t650t650) tNQUAD 
650  DO  66j  N»1,nPTS 
660  X (N)a2.*XPAX-X (N) 

CALL  line  (XfYtNPTSf ItOtO) 

DO  67o  N■l,^PTS 
670  Y (N)«2.»YPAx-Y (N) 

CALL  line  (XfY*NPISflt0»0) 

680  DO  69j  NxlfNPTS 


I 

i 


I 
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690  * (N)«^.*XKAX-X (N) 

CALL  line  (X,Y,NPIStl,0.0) 

700  CCNTInOE 

CALL  HLUT(l£.0*0*w«<»0) 

STOP 
END 

SLBPOuTINE  uRAm  (^tlL«JLtII«JjtXV(YV<Vf ISCAN) 

SCANS  GRID  and  plots  LINES  OF  EOOAL  VALOE 
DIMENblUN  Xv(l00)«TV(20U)tV(l00*^00) 

COMmOn/ORAMaRG/X (bOO) *7(6^0) tXMAX, YMAXyNQuAOtF IRSTX« GEL  TAX, FIRST Yt 
/DELTAy 
NT*0 
IXSO 
T»/ 

IA»IL 
00  3 
JAsJ 

AsV<lAf JA) 

f 0*Y(la»JA*I) 

IF  (A.EQ.C)  GO  TO  3 
IF  (A,LE.C..OR.D.LE.O.)  GOTO  3 
IF  (A.EO.T)  1,2 

1 T«A*(u-A)*I.E-6 

2 CALL  OA(OiA,T«XV,YVtX, Y,IA, JA,IX,NT, JJ,V) 

IF  (NT.nE.O)  go  to  30 

3 CCNTliNOE 
JA»JL 

00  6 l«ltll 

IA»IL-I 

Bsv ( lA« 1, JA) 

Asy ( I A, JA) 

IF  (B.EO.A)  go  to  6 
IF  (a.le,o..or.u.le,o.)  goto  6 
IF  (B«E0«T)  4»5 


4 T«H* (A-B)«1,E-6 

5 CALL  Aa(A»ri,TtXV,YV.X,Y,lA,jA,IX,NT,II,V) 
IF  (Nl  .NE*0)  GO  TO  30 

6 CONTINUE 
lAsO 

00  1C  Jai*JJ 
JA»JL-J 
C«V(1,JA*1) 
bsy ( 1 , JA) 

IF  (C.EO.B)  go  to  10 
IF  (Ja,E0.1)  goto  7 
IF  (C.LE.O..ORtB.LE.O.)  GOTO  10 

7 IF  (C.LE.O,.OR.B.LT.O.)  GOTO  10 
IF  (C.EO.T)  8t9 

8 TaC*<d-C)*l.E-6 

9 CALL  BC<8tC,T,XV,YV,X,Y,IA,jA,IX,NT,JJ,V) 
IF  (Nf.NE.o)  GO  To  30 

10  CONTINUE 
JA«0 

DO  lA  I»ltII 
IA»I 


I 


Csv(U*l«l) 

IF  (U.E^J,C)  60  TO  )<* 

IF  (Iu*EO>l.ANO*IbCAN.eU.?)  60TC  11 
IF  (C.Le.0.,OR.D.LE*O.)  GOTO  lA 
11  IF  (C.LF..U..OR.O.LT.O.  ) GUTO  14 
IF  (O.EO.T)  I2tl3 
1?  TsO*<C-i))*l.E-6 

13  CALL  CO(C»0,TtXV,YV,X,Y,lA,jA,lXtNT,II,V) 
IF  (NT.NE.O)  GO  To  30 

14  CCNTl^Ut 

IF  ( laCAN.iME.  1 ) GOTO  22 

DC  21  1=) »II 

DC  21  J*l»Jj 

IA»I 

JA*J 

C=v (l*i»j»i) 
b=v (1*1 t J) 

IF  (C.EQ.b)  GO  TO  le 
IF  (J.EO.l)  GOTO  15 
IF  (C.Lb.O>.OR.6.LF*0.)  GOTO  18 

15  IF  (C.Lt,C..OR.8.LT.O.)  GOfO  18 
IF  (C.EO.I)  16*17 

16  TsC*(d-C)«l,E-6 

17  CALL  dC(B.C,T,XV,YV»X,Y,lA,JA,IXtM,JJ,V) 
IF  (M.NF.U)  60  TO  3o 

18  D=V(I*J*1) 

IF  (C.EO.C)  60  TO  21 
IF  (C.Lt.O..OR.D.LE.O.)  GOTO  21 
IF  (C.EU.l)  1*5*20 

19  T=C*(u-C)*l.E-6 

20  CALL  CD«C*D,T.XV,YV*X,Y.lAtJAtIX.NT,IIfV) 
IF  (NT.MF.O)  GO  TO  30 

21  CONTINUE 
ISCANso 
GO  TO  69 

22  DC  29  J»l*Jj 
DO  29  I»1*I1 
IA  = I 

JA»J 

C=V<I*1*J»1) 

b=V(l*l*J) 

IF  (C.EU.B)  GO  TO  25 
IF  (C.LF.0t.OR.B.LE*O.)  GOTO  25 
IF  (C.EO.T)  23*24 

23  T=C*(8-C)»1.E-6 

24  CALL  dC(8*CtT*XVfYV*X,Y*lAtJA*IX*M,JJ»V) 
IF  (M.NE.O)  GO  To  30 

25  0«V(I*J*1> 

IF  (C.EO.O)  GO  TO  29 
IF  (l.EU.l)  GOTO  26 
IF  (C.LE.0*.OR.O«LE*0.)  GOTO  29 

26  IF  (C.LE.Ot.OR.O.LTiO.)  GOTO  29 
IF  (C.EQ.T)  27*28 

27  T»C*<U-C»*I.E-6 

28  CALL  CO(C*OiT*XV*YV*X»Y*IA,jA,lX*M,II*V) 
IF  (N| .NE*0)  GO  TO  30 
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29  CONTlNOt 
ISCAN«0 
6C  TO  69 

30  IST«U 

jst«J4 

31  IF  (U.EO.0.O9.IA.E(^>IL.O»4.JA.eG.0>OR«JA.£O.JL)  GO  TO  So 
NlXalX 

IF  (T.EO.V ( lAt JA) ) V(IA.JA)»V(IA.JA»*1.0001 

IF  (T.EU.V(IA«1«JA>)  V(lA«lf JA)3V(lA*lt jA)*1.0001 

IF  (T.EO.y (IA*1»JA*1) ) V(lA*i»JA*l)»V(IA*l.JA^l)*1.000l 

IF  (T.EU.y ( lAt JA*1 ) » V(IA,jA*l)»y(IAiJA»l)*l.0001 

Asy ( I Ai JA) 

bBy ( I A*i t JA) 

C«V (1 A*1 t JA*1) 

0»V  < Ixt JA*  1 ) 

IF  (lA.tO.l.ANO.UCAN.EU.a)  GOTO  32 
IF  (Ja.EQ.1)  goto  33 

IF  (A.lE.O>.OR.B.LE>O..OR.C*LE.O*.UR.O.LE.O<)  goto  40 

32  IF  (A.lT.u>.OR.B«LE>C..UR.C«I.E.O<.OR.O.LT.O.)  goto  so 
GOTO  34 

33  IF  (A.LT.O>.OR.H.LTfO<>OR.C>l>E.O*<OR.O.LE.O«)  GOTO  50 

34  IF  (Nl .EO.l)  GO  To  3b 

IF  (Nr.EO.2)  GO  TO  36 

IF  (Nl.EO.J)  GO  TO  37 

IF  (Nr.EO.4)  GO  TO  3U 

GO  TO  69 

35  Call  u a ( d » a * t t X y t t y t X I y * I a t ja f I X t n t * j j » V ) 

IF  (Ix.GT.MIX)  GOIO  39 

CALL  oC  (HtC,TtXy,Yy»X| Y«lA,JA,lXtNTt JJ«V) 

IF  (Ix.GT.NlX)  GOIO  39 

CALL  A0IA*B*TtXy*YV»XtY»lAtJA»IXtNT*II*y) 

GC  TO  39 

36  CALL  AS ( A t B f T f X y » y V » X t Y * I A * J A * I X * ^T * 1 1 1 y > 

IF  (Ix.GT.NlX)  GOIO  39 

CALL  CDlCtO.TtXytTy.X.Y.lA.JA.IX.NT.lI.V) 

IF  (U.GT.NlX)  GOIO  39 

CALL  3CIH»C.T.Xv,YytX,Y,IA,JA,IX*NT,JJ.y) 

GC  TO  39 

37  CALL  bC  l«tC,T.Xy,YV.X,Y,lA,JA,lX,NT,JJ,y) 

IF  (U.GT.miX)  goto  39 

CALL  LiA«ntA,T.Xy,YytX,Y*IA,jA,IXtNT.JJ*V) 

IF  (U.GT.NIX)  GOIO  39 

CALL  CD(C«0,T,Xy,YV.X,Y,lA,JA,lX,NT.IIiy) 

GC  TO  39 

3fl  CALL  CO<CtO,TiXy,YVtX,Y.lA,JA,IXfNT,II.y) 

IF  (lA.GT.NIX)  GOIO  39 

CALL  ASlAtS.TtXy.YytX.YflA.JA.IXtNTflltV) 

IF  (Ia.GT.nIX)  GOIO  39 

CALL  uAln*A,T»Xy»yv*X,Y,lA,JA,IX.NT,JJ.y) 

39  IF  (Ia.EQ.IsT.ANO.JA.EO.JST)  50.31 

40  IF  (X (IX) .EL.X (IX-1 ) ) GOTO  50 
S«IY(IX)-T(IX-1))/(X(IX)-X(IX-1)) 

SB*T(IX)-5*xIIX) 

41  SL» ITy I JA*1) -Yy ( JA) )/ (Xy ( IA)-xy ( 14»1) ) 

SR»(Yy ( jA*i) -yy ( JA) )/ (xy I iA*i) -xy  i ia) ) 
orthl»ass(Sl*i ./S) 

ORThRsASS ISr»1 ./S) 
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U (OrtTHH.LT  .OKThl)  (jOTO  42 
S9D»Tv(Ja)-SL*'‘V(U*1) 

XR» (baO-SB! / (S-SL) 

GOTO  43 

1.?  SBD»Tv(JA)-SH*X''tlA> 

XH*(SoD-SB)/ (S-SH» 

43  1XSIX*1 
X (I  X) =XH 

Y ( IX)  =S*X  ( IX)  *Stl 

50  CALL  xPLOlY (X.Y.^.IX) 

IF  (Xi«iAX.C-T.YPAX)  goto  60 

IF  (Nt.UA0.ec.2»ANU.  (2.*XMAX)  •GT.Yt<AX)  GOTO  feO 
DO  54  I»ltIX 

GCTO  (51.53,52*52) .NUlAD 

51  X(I)=XMAX-X(I) 

GCTO  34 

5?  Y (I)*YMAX-Y(I) 

33  X(I)=aMAX*X(I) 

34  CCNTI^UE 

X (IX*  I ) =F IHSTX 
X(lX*2l=0ELrAX 
Y(IX*l)aFIRSTY 
Y (IX*2)=DElTAY 
CALL  lINE(Y,X,IX,1,0,0) 

GCTo  (69,58,55,55) .NOUAD 
55  DO  56  l»l,lx 

36  X(I)=2.*XI'Ax-X(I) 

CALL  LII^E  (Y,X,IX,I  ,0,0) 

DO  57  1 = 1, IX 

37  Y(I)*2.*Y)'Ax-Y(I) 

CALL  lINE<Y,X, IX, 1,0,0) 

38  DC  59  1 = 1, Ix 

39  X(I)*2.*XI«4X-X(I) 

CALL  LII''E(Y,X, IX, 1,0,0) 

GCTo  69 

60  GCTo  (63,61,61 t61) ,NUUA0 

61  DO  62  I»1*IX 

X ( I ) =XMAX*X ( I ) 

62  IF  (NhUA0.£(;.4)  Y(I)=2,»YMaX-Y(I) 

63  X(IX*1)»FIHSTX 
X(1X*2)=OElTAX 
Y(IX*1)=FIRSTY 
Y(IX*2)»OElTAY 

CALL  line (X,Y,IX,1 ,0,0) 

GCTO  (69,67,64*64) ,NUUAO 

64  DO  65  I»1*IX 

65  X(I)*2.•X^'AX-X(I) 

CALL  LlNE(X,Y, IX, 1,0,0) 

DC  66  1 = 1, IX 

66  Y (1  )»<;.»Yr<AX-Y  (1) 

CALL  line  (X,Y,IX,1,(J,0) 

67  DO  68  1 = 1, IX 

68  X (I)*2,*XWAX-X  (I) 

CALL  lINEIX,Y,1X, 1,0,0) 

69  CCNTI6UE 

return 

END 
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SCBHUuTINE  RSOLV  (RiH4tHB*A2iAliA0) 

C 

C calculates  poots  uf  uuadhatic  equation 

c 

P=(-A1  - 5UPT(AI**2  - A.*A2*Ao>1  /(2.*A2) 
print  1 .RtPAtRBi A2tAl tAU 

1 format  (4H  P x»6E16.B) 

IF  (R.GE.Ra.AND.R.LE.RB)  goto  2 

P=(-Al  ♦ SUPT(A1**2  - 4.*A2«Ao))  /(2.*A2) 

PRINT  1 ,H,RA,RdiA2fAl,A0 

2 continue 

PETUMn 

END 

SLHROuTiNE  XPLOTY(X.Y,2,IX) 

C 

C PRiNli  COCRcINATES  OF  EOUaL  YaLLE  LINES 

c 

DIMENSION  X(l) *Y  ( n 
PRINT  3 
print  2.  IX 
print  At  i 

print  It  (X  (K)  t N=1 1 IX) 

PRiNt  3 

PRINT  It  (Y(K).  N=ltIX) 
print  3 

1 FCRMATdH  ,dEl6.8) 

FCPMATdH  tAl6t7H  POINTS) 

FCRMATdH  ) 

FCHmAI(9h  value  = tF8,3) 

return 

END 

Subroutine  AB(AtBtTtxv*YVtXtY»iAtjAtiXtNT*iitV) 

C searches  dOTTOM  side  CF  rectangle  for  given  value 

c 

dimension  Xvd)tYvd)tXd)tY(l)tVdOOt20  0> 

IF  (A.GT.T.ANDtB.LT.T.OR.A.LT.T.AND.B.GT.T)  it? 

1 Pl»XVdA) 

01  = XV ( IA*1) 

XX=<PltQl)/2. 

T1=A 

TEsR 

IF  (lA.EOtl)  2i3 

2 Rl  = XVdA*2) 

T3*V ( I A*?t JA) 

CALL  SOL33(XX,TtPltul,RltTltT2,T3) 

GC  TO  6 

3 IF  (Ia.EO.II)  At5 

<*  P)=XV(IA-1) 

T3»V (lA-l t JA) 

CALL  SOL33(XX,T,PltultRltTltl2tT3) 

GC  TO  6 

S IF  (V lIA-l,jA) .LEtO? I GOTO  2 
IF  (V  dA*2t  jA)  .LE.Ct)  GOTO  A 
Rl  = XVdA-l) 

5l=XV ( IA*2) 

T3=V(IA-ltJA) 
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T4*V(iA*2tJ/>) 

CALL  sOLAA  (XX,T  tPlfUl  *S1  tTl»T2»T3«TA) 
b NT31 

IX»lX*i 
X ( IX) SXX 
Y ( IX)=YV ( JA) 

JA3JA-1 
7 PFTuRn 
EKO 

SLBROuTiNE  BC(B.C.T»Xv»YV,XtY,IAtjA,lX,NT,JJtV) 
C 

C SEARCHES  R.H,  Slot  OF  RECTaNOLE  FOR  GIVEN  VALUE 
C 

dimension  Xv(l),YV(l),X(l),Y(l),V(100t200) 

IF  <B,GI.T.aN0.C.LT.T.OR,h,lI,T.aM).C.GT.T)  li7 
1 P1=YV(JA) 

U1=YV ( JA*1) 

XXs(Pi*Jl)/2. 

T1  = H 
T2  = C 

IF  (JA.EQ.I)  2«3 
^ Rl=YV(JA*2) 

T3=v<iA»l *JA*2) 

CALL  SOL33 (XX,T*P1 lUl f RltTl . I2«T3) 

GC  TO  6 

i IF  A»5 

A R1*YV(JA-1) 

T3*V(IA+ltJA-l) 

CALL  SOL33(XX,T,Pl»Ul,R),TltT2tT3) 

GC  TO  b 

b IF  (V  (lA*l,jA-l) .LE«0*)  GOTO  2 
IF  (V(IA*ltjA*2) .LE.O.)  GOTO  4 
Rl*YV ( jA-1) 

S13YV { JA*2) 

T3aV(IA*l»jA-l) 

T4=V(1A*1*JA*2) 

CALL  bOL44 (XX,T,Pl,ai,Rl,Slf Tj,T2,T3iT4) 
b NT^a 
IX=IX*1 

X ( IX) =XV ( I4»l ) 

Y (1X)=XX 
lAsIA*! 

7 return 
END 

SLBROyTlNE  CO(C*0»T*Xy,YV,X,Y,IA«JA,IX,NT,II»V) 

searches  Top  side  of  rectangle  for  given  value 

DIMENSION  XV  n ) yYV (1) *X(I) »Y (1) tV (100*200) 

IF  (C.6T.T.aND>D*LT.T.0R.C.LT.T.AN0.0>6T.T>  1*7 

1 Pl^XVdA) 

(J1*XV(IA*1) 

XX»(Pl*Ul)/2. 

T1*D 
T2=C 

IF  (Ia.EO.1)  2*3 

2 R1>xV(IA«2) 
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T3»V ( lA*2t JA*1 > 

CALL  SOL33  (xX,T,Pi,(Jl,ai,TliT2,TJ) 
ec  TO  6 

3 IF  (iM.bO.Il)  A«5 

4 Pl«XV(IA-l) 

T3»V (lA-1 * JA41 ) 

CALL  bOL33(xX,T,PlfUl,aitTl»T2»T3) 

GC  TO  6 

b IF  .LE.O.)  GOTO  2 

IF  (>/(IA*2ijA*l)  .LE«0.)  GOTO  4 
Pl«XVaA-l) 

Sl*XV ( IA»2) 

T3sv(lA-ltJA4l) 

T4*V ( i A42f JA*1 ) 

CALL  bOL44  (XX,T,Pl .01 .«! *Sl *ll*T2.T3tT4» 
t NT  = 3 
IX=IX*1 
X(IX)*XA 
Y(IX)=YV<JA*1) 

JA»JA*1 
V PETOHN 
END 

SLBHOuTINE  L/A  (DtA.T»XV.YV,X*Y  tIAt  JA.IX.NT,  JJtV) 
SEAPCHES  L.H.  side  of  HECTaNGLE  FOP  GlvEN  VALUE 


dimension  Xv(1)*YV(1).X(1»,Y(1>*V(100»200) 

IF  (D.GT .T . aNO.A.LT.T.OP.D.LT .T.AND.A.GT .T)  1*7 

1 PISYV(JA) 

01=YVlJA*l) 

XX=  (P  l*Ol ) /2. 

T)*A 

T2=n 

IF  (Ja.EQ.1)  2*3 

2 P1«YV(JA*2) 

T3=V{iA*JA*2) 

CALL  SOL 33  lXX,TfPl*01,Rl»Tl*l2*T3) 

GC  TO  6 

3 IF  (Ja.EQ.Jj)  4.5 

4 Pl*YV(JA-l) 

T3SV ( lA. jA-l) 

CALL  SOL33(XX,T,Pl  ,01,Rl,Ti,T2,T3) 

GC  TO  6 

b IF  (V  (IA,JA-1) .LE.O*)  GOTO  2 
IF  (V  (IA,JA*2) .LE.O*)  GOTO  4 

PiaYV(J4-l) 

Sl*YV ( JA*2) 

T3*V ( 1 A. jA-\ ) 

T4»V (IA*JA*2) 

CALL  SOL44(XX.T.Pl*01.Rl*Sl*ll*T2*T3.T4> 
e NTa4 
IXsIX.l 
X(lX)sXV(lA) 

Y(1X)=XX 
IA»IA-1 
7 PETURfN 
END 
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StHHUuflNt  SOL33 (AXi T*P1 i0lt«l*TltT2.T3> 

polynomul  interpolation  »»ITM  3 ncoes 

XXsXX-Pl 

XISQ. 

X2sui-Hl 

X3»R1-P1 

PQ  = X2 

XX2»X£'*X2 

XX3»XJ*X3 

TT*T-  I I 

TT2»T2-T 1 

TT3*T3-I1 

D7»XX2*X3-Xa3»X2 

AA»(T12*X3-IT3*<2)/UT 

bB= (Xx?*TT3-XX3*T I 2) /OT 

TXs-Ti 

Ti=TTc-(T 

1 TY*(AA*Xx*cld)*XX-n 

tl®tx*ty 

IF  (Tu>  2*5*3 

2 X2=XX 
TZ«TY 
(jC  TO  4 

3 XlsxX 

TX3TT 

4 XLi=(Xi*X2)/2, 

0 = AijS(  IXX-Xl)/PO> 

XXSXO 

IF  (O,lI.*005)  b*l 

j,  XX*XX*P1 

petoRn 

EFO 

SLHROoTInE  S0L44 (xX,1 ,P1,Q1 ,Ri,S1 *Ti *I2*T3*T4) 
RCLYNyMlAL  INTERPOLATION  **ITH  4 nOOES 
XXaXX-Pl 

X1*0* 

X2S(JI-P1 

X3*R1-Pl 

X4»Sl-Pl 

PG*X2 

XX2»Xi;*X? 

XX3aXj*X3 

XX4aX4«X4 

XXX2SAX24X2 

XXX3*xX34X3 

XXX43XX4*X4 

TT«T-f  1 
TT2»T2-T1 
TT3«t3-|l 
T7  4»f4**f  1 

DT*Dt3(XXX2*XX2*Xt:tXXX3*XX3*X3,XXX4,XX4*X4) 
AA»DT3 ( IT2*XX2*X2*TT3*XX3,x3*TTA*XX4*X4) /DT 
BB»OTj(XXX2.TT2*X2*XXX3*TT3»X3*XXX4,TT4*X4)/OI 


COPY 


cc*0f  J I2*;ixX3.XX3»TT3»KXX4»XA4.TT4)/L)T 

Ixs-Tl 

T2  = TTe-n 

1 TV=( iAA*XX*HB)*XX*CC)*XX-TT 
TOsTX*TY 
IF  (To)  ?*5t3 
^ XjsXX 
T2»TY 
(iC  TO  4 

3 XJSXX 
TX  = TV 

4 Xt*(Xl»X?)/^. 

0=ABi>(  (XX-xo)  /MU) 

XX*XU 

IF  (U,Ul.*i)05)  5*1 

t XX*XX*pi 

ENO 

function  Cl3 (A1 *Bl *Cl,A2,B2*C2,A3.B3iC3) 

C 

C CALCOlAIES  uETEHMINANT  of  3 By  3 matrix 

c 

DT3aAi*b?*C3*A2*B3*Cl*A3*Hi*C2-Al*83*C2-A2*Bi*C3-A3*B2*Cl 

return 

END 
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